From efadc84593ba47aa944e1caef6d3b8aa70cdd6c8 Mon Sep 17 00:00:00 2001
From: TimQu <tim.quatmann@rwth-aachen.de>
Date: Mon, 10 Aug 2015 19:56:14 +0200
Subject: [PATCH] Beautified the Code, removed unused stuff, minor improvements

Former-commit-id: 4c16f9163c8d639a88bc87673ecda9589cfeefb2
---
 src/adapters/CarlAdapter.h                    |    4 +-
 src/adapters/Smt2ExpressionAdapter.h          |    4 +-
 .../SparseDtmcRegionModelChecker.cpp          | 1605 -----------------
 .../SparseDtmcRegionModelChecker.h            |  386 ----
 .../region/ApproximationModel.cpp             |   14 +-
 src/modelchecker/region/ApproximationModel.h  |    3 +-
 src/modelchecker/region/ParameterRegion.cpp   |  258 +++
 src/modelchecker/region/ParameterRegion.h     |  136 ++
 src/modelchecker/region/SamplingModel.cpp     |    8 +-
 src/modelchecker/region/SamplingModel.h       |    3 +-
 .../region/SparseDtmcRegionModelChecker.cpp   |  734 ++++++++
 .../region/SparseDtmcRegionModelChecker.h     |  205 +++
 src/solver/Smt2SmtSolver.cpp                  |    6 +-
 src/utility/cli.h                             |   13 +-
 src/utility/regions.cpp                       |  155 +-
 src/utility/regions.h                         |  123 +-
 16 files changed, 1428 insertions(+), 2229 deletions(-)
 delete mode 100644 src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp
 delete mode 100644 src/modelchecker/reachability/SparseDtmcRegionModelChecker.h
 create mode 100644 src/modelchecker/region/ParameterRegion.cpp
 create mode 100644 src/modelchecker/region/ParameterRegion.h
 create mode 100644 src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
 create mode 100644 src/modelchecker/region/SparseDtmcRegionModelChecker.h

diff --git a/src/adapters/CarlAdapter.h b/src/adapters/CarlAdapter.h
index 3100b13f8..14bfaaf07 100644
--- a/src/adapters/CarlAdapter.h
+++ b/src/adapters/CarlAdapter.h
@@ -10,7 +10,7 @@
 #include <carl/core/MultivariatePolynomial.h>
 #include <carl/core/RationalFunction.h>
 #include <carl/core/VariablePool.h>
-#include <carl/core/Constraint.h>
+#include <carl/formula/Constraint.h>
 #include <carl/core/FactorizedPolynomial.h>
 
 namespace carl {
@@ -39,7 +39,7 @@ namespace storm {
         typedef cln::cl_RA Coefficient;
     typedef carl::MultivariatePolynomial<Coefficient> RawPolynomial;
     typedef carl::FactorizedPolynomial<RawPolynomial> Polynomial;
-	typedef carl::CompareRelation CompareRelation;
+	typedef carl::Relation CompareRelation;
 	typedef carl::RationalFunction<Polynomial> RationalFunction;
 }
 
diff --git a/src/adapters/Smt2ExpressionAdapter.h b/src/adapters/Smt2ExpressionAdapter.h
index f0d9cea08..bc6c0807f 100644
--- a/src/adapters/Smt2ExpressionAdapter.h
+++ b/src/adapters/Smt2ExpressionAdapter.h
@@ -68,7 +68,7 @@ namespace storm {
              */
             std::string translateExpression(carl::Constraint<storm::RawPolynomial> const& constraint) {
                                 
-                return  "(" + carl::toString(constraint.rel()) + " " +
+                return  "(" + carl::relationToString(constraint.relation()) + " " +
                             constraint.lhs().toString(false, useReadableVarNames) + " " +
                             "0 " +
                         ")";
@@ -81,7 +81,7 @@ namespace storm {
              */
             std::string translateExpression(carl::Constraint<storm::Polynomial> const& constraint) {
                                 
-                return  "(" + carl::toString(constraint.rel()) + " " +
+                return  "(" + carl::relationToString(constraint.relation()) + " " +
                             constraint.lhs().toString(false, useReadableVarNames) + " " +
                             "0 " +
                         ")";
diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp
deleted file mode 100644
index 66378024b..000000000
--- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp
+++ /dev/null
@@ -1,1605 +0,0 @@
-#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
-
-#include <chrono>
-
-#include "src/adapters/CarlAdapter.h"
-
-
-#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
-#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
-
-#include "src/utility/graph.h"
-#include "src/utility/vector.h"
-#include "src/utility/macros.h"
-
-#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
-#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-
-#include "src/modelchecker/region/ApproximationModel.h"
-#include "src/modelchecker/region/SamplingModel.h"
-
-#include "src/exceptions/InvalidPropertyException.h"
-#include "src/exceptions/InvalidStateException.h"
-#include "src/exceptions/UnexpectedException.h"
-
-
-namespace storm {
-    namespace modelchecker {
-        
-        
-        template<typename ParametricType, typename ConstantType>
-        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::ParameterRegion(std::map<VariableType, CoefficientType> lowerBounds, std::map<VariableType, CoefficientType> upperBounds) : lowerBounds(lowerBounds), upperBounds(upperBounds), checkResult(RegionCheckResult::UNKNOWN) {
-            //check whether both mappings map the same variables and precompute the set of variables
-            for(auto const& variableWithBound : lowerBounds) {
-                STORM_LOG_THROW((upperBounds.find(variableWithBound.first)!=upperBounds.end()), storm::exceptions::InvalidArgumentException, "Couldn't create region. No upper bound specified for Variable " << variableWithBound.first);
-                this->variables.insert(variableWithBound.first);
-            }
-            for(auto const& variableWithBound : upperBounds){
-                STORM_LOG_THROW((this->variables.find(variableWithBound.first)!=this->variables.end()), storm::exceptions::InvalidArgumentException, "Couldn't create region. No lower bound specified for Variable " << variableWithBound.first);
-            }
-            
-        }
-
-        template<typename ParametricType, typename ConstantType>        
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setViolatedPoint(std::map<VariableType, CoefficientType> const& violatedPoint) {
-            this->violatedPoint = violatedPoint;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getViolatedPoint() const {
-            return violatedPoint;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setSatPoint(std::map<VariableType, CoefficientType> const& satPoint) {
-            this->satPoint = satPoint;
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getSatPoint() const {
-            return satPoint;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setCheckResult(RegionCheckResult checkResult) {
-            //a few sanity checks
-            STORM_LOG_THROW((this->checkResult==RegionCheckResult::UNKNOWN || checkResult!=RegionCheckResult::UNKNOWN),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from something known to UNKNOWN ");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSSAT || checkResult!=RegionCheckResult::EXISTSVIOLATED),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSSAT to EXISTSVIOLATED");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSSAT || checkResult!=RegionCheckResult::ALLVIOLATED),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSSAT to ALLVIOLATED");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSVIOLATED || checkResult!=RegionCheckResult::EXISTSSAT),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSVIOLATED to EXISTSSAT");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSVIOLATED || checkResult!=RegionCheckResult::ALLSAT),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSVIOLATED to ALLSAT");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSBOTH || checkResult!=RegionCheckResult::ALLSAT),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSBOTH to ALLSAT");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::EXISTSBOTH || checkResult!=RegionCheckResult::ALLVIOLATED),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSBOTH to ALLVIOLATED");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::ALLSAT || checkResult==RegionCheckResult::ALLSAT),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from ALLSAT to something else");
-            STORM_LOG_THROW((this->checkResult!=RegionCheckResult::ALLVIOLATED || checkResult==RegionCheckResult::ALLVIOLATED),storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from ALLVIOLATED to something else");
-            this->checkResult = checkResult;
-        }
-
-        template<typename ParametricType, typename ConstantType>        
-        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::RegionCheckResult SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getCheckResult() const {
-            return checkResult;
-        }
-
-      
-                
-        template<typename ParametricType, typename ConstantType>
-        std::set<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getVariables() const{
-            return this->variables;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getLowerBound(VariableType const& variable) const{
-            auto const& result = lowerBounds.find(variable);
-            STORM_LOG_THROW(result!=lowerBounds.end(), storm::exceptions::IllegalArgumentException, "tried to find a lower bound of a variable that is not specified by this region");
-            return (*result).second;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getUpperBound(VariableType const& variable) const{
-            auto const& result = upperBounds.find(variable);
-            STORM_LOG_THROW(result!=upperBounds.end(), storm::exceptions::IllegalArgumentException, "tried to find an upper bound of a variable that is not specified by this region");
-            return (*result).second;
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        const std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getUpperBounds() const {
-            return upperBounds;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        const std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getLowerBounds() const {
-            return lowerBounds;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        std::vector<std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType>> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getVerticesOfRegion(std::set<VariableType> const& consideredVariables) const{
-            std::size_t const numOfVariables=consideredVariables.size();
-            std::size_t const numOfVertices=std::pow(2,numOfVariables);
-            std::vector<std::map<VariableType, CoefficientType>> resultingVector(numOfVertices,std::map<VariableType, CoefficientType>());
-            if(numOfVertices==1){
-                //no variables are given, the returned vector should still contain an empty map
-                return resultingVector;
-            }
-            
-            for(uint_fast64_t vertexId=0; vertexId<numOfVertices; ++vertexId){
-                //interprete vertexId as a bit sequence
-                //the consideredVariables.size() least significant bits of vertex will always represent the next vertex
-                //(00...0 = lower bounds for all variables, 11...1 = upper bounds for all variables)
-                std::size_t variableIndex=0;
-                for(auto const& variable : consideredVariables){
-                    if( (vertexId>>variableIndex)%2==0  ){
-                        resultingVector[vertexId].insert(std::pair<VariableType, CoefficientType>(variable, getLowerBound(variable)));
-                    }
-                    else{
-                        resultingVector[vertexId].insert(std::pair<VariableType, CoefficientType>(variable, getUpperBound(variable)));
-                    }
-                    ++variableIndex;
-                }
-            }
-            return resultingVector;
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        std::string SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::checkResultToString() const {
-            switch (this->checkResult) {
-                case RegionCheckResult::UNKNOWN:
-                    return "unknown";
-                case RegionCheckResult::EXISTSSAT:
-                    return "ExistsSat";
-                case RegionCheckResult::EXISTSVIOLATED:
-                    return "ExistsViolated";
-                case RegionCheckResult::EXISTSBOTH:
-                    return "ExistsBoth";
-                case RegionCheckResult::ALLSAT:
-                    return "allSat";
-                case RegionCheckResult::ALLVIOLATED:
-                    return "allViolated";
-            }
-            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not identify check result")
-            return "ERROR";
-        }
-            
-        template<typename ParametricType, typename ConstantType>
-        std::string SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::toString() const {
-            std::stringstream regionstringstream;
-            for(auto var : this->getVariables()){
-                regionstringstream << storm::utility::regions::convertNumber<SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType,double>(this->getLowerBound(var));
-                regionstringstream << "<=";
-                regionstringstream << storm::utility::regions::getVariableName(var);
-                regionstringstream << "<=";
-                regionstringstream << storm::utility::regions::convertNumber<SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType,double>(this->getUpperBound(var));
-                regionstringstream << ",";
-            }
-            std::string regionstring = regionstringstream.str();
-            //the last comma should actually be a semicolon
-            regionstring = regionstring.substr(0,regionstring.length()-1) + ";";
-            return regionstring;
-        }
-                
-        
-        
-        template<typename ParametricType, typename ConstantType>
-        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc<ParametricType> const& model) : 
-                model(model),
-                eliminationModelChecker(model),
-                smtSolver(nullptr),
-                probabilityOperatorFormula(nullptr),
-                samplingModel(nullptr),
-                approximationModel(nullptr),
-                isReachProbFunctionComputed(false),
-                isResultConstant(false) {
-            //intentionally left empty
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::canHandle(storm::logic::Formula const& formula) const {
-             //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ])
-            if(!formula.isStateFormula()){
-                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected a stateFormula");
-                return false;
-            }
-            if(!formula.asStateFormula().isProbabilityOperatorFormula()){
-                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected a probabilityOperatorFormula");
-                return false;
-            }
-            storm::logic::ProbabilityOperatorFormula const& probOpForm=formula.asStateFormula().asProbabilityOperatorFormula();
-            if(!probOpForm.hasBound()){
-                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: The formula has no bound");
-                return false;
-            }
-            if(!probOpForm.getSubformula().asPathFormula().isEventuallyFormula()){
-                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected an eventually subformula");
-                return false;
-            }
-            if(model.getInitialStates().getNumberOfSetBits() != 1){
-                STORM_LOG_DEBUG("Input model is required to have exactly one initial state.");
-                return false;
-            }
-            return true;
-        }
-            
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::specifyFormula(storm::logic::Formula const& formula) {
-            std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now();
-            STORM_LOG_THROW(this->canHandle(formula), storm::exceptions::IllegalArgumentException, "Tried to specify a formula that can not be handled.");
-            
-            //Get subformula, target states
-            //Note: canHandle already ensures that the formula has the right shape and that the model has a single initial state.
-            this->probabilityOperatorFormula= std::unique_ptr<storm::logic::ProbabilityOperatorFormula>(new storm::logic::ProbabilityOperatorFormula(formula.asStateFormula().asProbabilityOperatorFormula()));
-            storm::logic::EventuallyFormula const& eventuallyFormula = this->probabilityOperatorFormula->getSubformula().asPathFormula().asEventuallyFormula();
-            std::unique_ptr<CheckResult> targetStatesResultPtr = this->eliminationModelChecker.check(eventuallyFormula.getSubformula());
-            storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector();
-            
-            computeSimplifiedModel(targetStates);
-            initializeSampleAndApproxModel();
-           
-            //some information for statistics...
-            std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now();
-            this->timePreprocessing= timePreprocessingEnd - timePreprocessingStart;
-            this->numOfCheckedRegions=0;
-            this->numOfRegionsSolvedThroughSampling=0;
-            this->numOfRegionsSolvedThroughApproximation=0;
-            this->numOfRegionsSolvedThroughSubsystemSmt=0;
-            this->numOfRegionsSolvedThroughFullSmt=0;
-            this->numOfRegionsExistsBoth=0;
-            this->numOfRegionsAllSat=0;
-            this->numOfRegionsAllViolated=0;
-            this->timeCheckRegion=std::chrono::high_resolution_clock::duration::zero();
-            this->timeSampling=std::chrono::high_resolution_clock::duration::zero();
-            this->timeApproximation=std::chrono::high_resolution_clock::duration::zero();
-            this->timeMDPBuild=std::chrono::high_resolution_clock::duration::zero();
-            this->timeSubsystemSmt=std::chrono::high_resolution_clock::duration::zero();
-            this->timeFullSmt=std::chrono::high_resolution_clock::duration::zero();
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::computeSimplifiedModel(storm::storage::BitVector const& targetStates){
-                        
-            //Compute the subset of states that have a probability of 0 or 1, respectively and reduce the considered states accordingly.
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->model, storm::storage::BitVector(this->model.getNumberOfStates(),true), targetStates);
-            storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
-            storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
-            storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
-            // If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction.
-            if (this->model.getInitialStates().isDisjointFrom(maybeStates)) {
-                STORM_LOG_WARN("The probability of the initial state is constant (0 or 1)");
-                this->reachProbFunction = statesWithProbability0.get(*(this->model.getInitialStates().begin())) ? storm::utility::zero<ParametricType>() : storm::utility::one<ParametricType>();
-                this->isReachProbFunctionComputed=true;
-                this->isResultConstant=true;
-            }
-            // Determine the set of states that is reachable from the initial state without jumping over a target state.
-            storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->model.getTransitionMatrix(), this->model.getInitialStates(), maybeStates, statesWithProbability1);
-            // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
-            maybeStates &= reachableStates;
-            // Create a vector for the probabilities to go to a state with probability 1 in one step.
-            std::vector<ParametricType> oneStepProbabilities = this->model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
-            // Determine the initial state of the sub-model.
-            storm::storage::sparse::state_type initialState = *(this->model.getInitialStates() % maybeStates).begin();
-            // We then build the submatrix that only has the transitions of the maybe states.
-            storm::storage::SparseMatrix<ParametricType> submatrix = this->model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
-            
-            // eliminate all states with only constant outgoing transitions
-            //TODO: maybe also states with constant incoming tranistions. THEN the ordering of the eliminated states does matter.
-            std::chrono::high_resolution_clock::time_point timeInitialStateEliminationStart = std::chrono::high_resolution_clock::now();
-            // Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
-            FlexibleMatrix flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix);
-            FlexibleMatrix flexibleBackwardTransitions= this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix.transpose(), true);
-            // Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated.
-            storm::storage::BitVector subsystem = storm::storage::BitVector(submatrix.getRowCount(), true);
-            //temporarily unselect the initial state to skip it.
-            subsystem.set(initialState, false);
-            this->hasOnlyLinearFunctions=true;
-            boost::optional<std::vector<ParametricType>> missingStateRewards;
-            for (auto const& state : subsystem) {
-                bool onlyConstantOutgoingTransitions=true;
-                for(auto const& entry : submatrix.getRow(state)){
-                    if(!this->parametricTypeComparator.isConstant(entry.getValue())){
-                        onlyConstantOutgoingTransitions=false;
-                        this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(entry.getValue());
-                    }
-                }
-                if(!this->parametricTypeComparator.isConstant(oneStepProbabilities[state])){
-                    onlyConstantOutgoingTransitions=false;
-                    this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(oneStepProbabilities[state]);
-                }
-                if(onlyConstantOutgoingTransitions){
-                    this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
-                    subsystem.set(state,false);
-                }
-            }
-            subsystem.set(initialState, true);
-            std::chrono::high_resolution_clock::time_point timeInitialStateEliminationEnd = std::chrono::high_resolution_clock::now();
-            this->timeInitialStateElimination = timeInitialStateEliminationEnd-timeInitialStateEliminationStart;
-            STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl);
-            std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl;
-            
-            //Build the simplified model
-            //The matrix. The flexibleTransitions matrix might have empty rows where states have been eliminated.
-            //The new matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
-            std::vector<storm::storage::sparse::state_type> newStateIndexMap(flexibleTransitions.getNumberOfRows(), flexibleTransitions.getNumberOfRows()); //initialize with some illegal index to easily check if a transition leads to an unselected state
-            storm::storage::sparse::state_type newStateIndex=0;
-            for(auto const& state : subsystem){
-                newStateIndexMap[state]=newStateIndex;
-                ++newStateIndex;
-            }
-            //We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead
-            storm::storage::sparse::state_type numStates=newStateIndex+2;
-            storm::storage::sparse::state_type targetState=numStates-2;
-            storm::storage::sparse::state_type sinkState= numStates-1;
-            //We can now fill in the data.
-            storm::storage::SparseMatrixBuilder<ParametricType> matrixBuilder(numStates, numStates);
-            for(storm::storage::sparse::state_type oldStateIndex : subsystem){ 
-                ParametricType valueToSinkState=storm::utility::regions::getNewFunction<ParametricType, CoefficientType>(storm::utility::one<CoefficientType>());
-                //go through columns:
-                for(auto const& entry: flexibleTransitions.getRow(oldStateIndex)){ 
-                    STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=flexibleTransitions.getNumberOfRows(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated.");
-                    valueToSinkState-=entry.getValue();
-                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],entry.getValue());
-                }
-                //transition to target state
-                if(!this->parametricTypeComparator.isZero(oneStepProbabilities[oldStateIndex])){ 
-                    valueToSinkState-=oneStepProbabilities[oldStateIndex];
-                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, oneStepProbabilities[oldStateIndex]);
-                }
-                //transition to sink state
-                if(!this->parametricTypeComparator.isZero(valueToSinkState)){ 
-                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, valueToSinkState);
-                }
-            }
-            //add self loops on the additional states (i.e., target and sink)
-            matrixBuilder.addNextValue(targetState, targetState, storm::utility::one<ParametricType>());
-            matrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one<ParametricType>());
-            //The labeling
-            storm::models::sparse::StateLabeling labeling(numStates);
-            storm::storage::BitVector initLabel(numStates, false);
-            initLabel.set(newStateIndexMap[initialState], true);
-            labeling.addLabel("init", std::move(initLabel));
-            storm::storage::BitVector targetLabel(numStates, false);
-            targetLabel.set(targetState, true);
-            labeling.addLabel("target", std::move(targetLabel));
-            storm::storage::BitVector sinkLabel(numStates, false);
-            sinkLabel.set(sinkState, true);
-            labeling.addLabel("sink", std::move(sinkLabel));
-            // other ingredients
-            boost::optional<std::vector<ParametricType>> noStateRewards;
-            boost::optional<storm::storage::SparseMatrix<ParametricType>> noTransitionRewards;  
-            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;            
-            // the final model
-            this->simplifiedModel = std::make_shared<storm::models::sparse::Dtmc<ParametricType>>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSampleAndApproxModel(){
-            
-
-            
-            //now create the models
-            this->approximationModel=std::make_shared<ApproximationModel>(*this->simplifiedModel);
-            this->samplingModel=std::make_shared<SamplingModel>(*this->simplifiedModel);
-            
-        }
-       /*
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSampleDtmcAndApproxMdp2(
-                    std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>>& sampleDtmc,
-                    std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>> & sampleDtmcMapping,
-                    std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> & approxMdp,
-                    std::vector<std::tuple<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&, std::size_t>> & approxMdpMapping,
-                    std::vector<std::map<VariableType, TypeOfBound>> & approxMdpSubstitutions,
-                    storm::storage::BitVector const& subsys,
-                    storm::storage::SparseMatrix<ParametricType> const& transitions,
-                    std::vector<ParametricType> const& oneStepProbs,
-                    storm::storage::sparse::state_type const& initState) {
-            
-            
-            //the matrix "transitions" might have empty rows where states have been eliminated.
-            //The DTMC/ MDP matrices should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
-            std::vector<storm::storage::sparse::state_type> newStateIndexMap(transitions.getRowCount(), transitions.getRowCount()); //initialize with some illegal index to easily check if a transition leads to an unselected state
-            storm::storage::sparse::state_type newStateIndex=0;
-            for(auto const& state : subsys){
-                newStateIndexMap[state]=newStateIndex;
-                ++newStateIndex;
-            }
-            //We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead
-            storm::storage::sparse::state_type numStates=newStateIndex+2;
-            storm::storage::sparse::state_type targetState=numStates-2;
-            storm::storage::sparse::state_type sinkState= numStates-1;
-            
-            //We can now fill in the dummy data of the transition matrices for dtmc and mdp
-            std::vector<ParametricType> oneStepToSink(oneStepProbs.size(), storm::utility::zero<ParametricType>()); // -2 because not needed for target and sink state
-            storm::storage::SparseMatrixBuilder<ConstantType> dtmcMatrixBuilder(numStates, numStates);
-            storm::storage::SparseMatrixBuilder<ConstantType> mdpMatrixBuilder(0, numStates, 0, true, true, numStates);
-            std::vector<std::map<VariableType, TypeOfBound> > mdpRowSubstitutions;
-            storm::storage::sparse::state_type currentMdpRow=0;
-            //go through rows:
-            for(storm::storage::sparse::state_type oldStateIndex : subsys){ 
-                ParametricType valueToSinkState=storm::utility::regions::getNewFunction<ParametricType, CoefficientType>(storm::utility::one<CoefficientType>());
-                // the dtmc and the mdp rows need to sum up to one, therefore the first entry that we add has value one.
-                ConstantType dummyEntry=storm::utility::one<ConstantType>(); 
-                // store the columns and values that we have added because we need the same information for the next rows in the mdp
-                std::vector<std::pair<storm::storage::sparse::state_type, ConstantType>> rowInformation;
-                mdpMatrixBuilder.newRowGroup(currentMdpRow);
-                std::set<VariableType> occurringParameters;
-                //go through columns:
-                for(auto const& entry: transitions.getRow(oldStateIndex)){ 
-                    valueToSinkState-=entry.getValue();
-                    storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringParameters);
-                    dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],dummyEntry);
-                    mdpMatrixBuilder.addNextValue(currentMdpRow, newStateIndexMap[entry.getColumn()], dummyEntry);
-                    rowInformation.emplace_back(std::make_pair(newStateIndexMap[entry.getColumn()], dummyEntry));
-                    dummyEntry=storm::utility::zero<ConstantType>(); 
-                    STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=transitions.getRowCount(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated.");
-                }
-                //transition to target state
-                if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ 
-                    valueToSinkState-=oneStepProbs[oldStateIndex];
-                    storm::utility::regions::gatherOccurringVariables(oneStepProbs[oldStateIndex], occurringParameters);
-                    dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, dummyEntry);
-                    mdpMatrixBuilder.addNextValue(currentMdpRow, targetState, dummyEntry);
-                    rowInformation.emplace_back(std::make_pair(targetState, dummyEntry));
-                    dummyEntry=storm::utility::zero<ConstantType>(); 
-                }
-                //transition to sink state
-                if(!this->parametricTypeComparator.isZero(valueToSinkState)){ 
-                    dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, dummyEntry);
-                    mdpMatrixBuilder.addNextValue(currentMdpRow, sinkState, dummyEntry);
-                    rowInformation.emplace_back(std::make_pair(sinkState, dummyEntry));
-                    oneStepToSink[oldStateIndex]=valueToSinkState;
-                }
-                // Find the different combinations of occuring variables and lower/upper bounds (mathematically speaking we want to have the set 'occuringVariables times {u,l}' )
-                std::size_t numOfSubstitutions=std::pow(2,occurringParameters.size()); //1 substitution = 1 combination of lower and upper bounds
-                for(uint_fast64_t substitutionId=0; substitutionId<numOfSubstitutions; ++substitutionId){
-                    //interprete substitutionId as a bit sequence
-                    //the occurringParameters.size() least significant bits of substitutionId will always represent the next substitution that we have to add
-                    //(00...0 = lower bounds for all parameters, 11...1 = upper bounds for all parameters)
-                    mdpRowSubstitutions.emplace_back();
-                    std::size_t parameterIndex=0;
-                    for(auto const& parameter : occurringParameters){
-                        if( (substitutionId>>parameterIndex)%2==0  ){
-                            mdpRowSubstitutions.back().insert(std::make_pair(parameter, TypeOfBound::LOWER));
-                        }
-                        else{
-                            mdpRowSubstitutions.back().insert(std::make_pair(parameter, TypeOfBound::UPPER));
-                        }
-                        ++parameterIndex;
-                    }
-                }
-                //We already added one row in the MDP. Now add the remaining 2^(occuringParameters.size())-1 rows
-                for(++currentMdpRow; currentMdpRow<mdpRowSubstitutions.size(); ++currentMdpRow){
-                    for(auto const& entryOfThisRow : rowInformation){
-                        mdpMatrixBuilder.addNextValue(currentMdpRow, entryOfThisRow.first, entryOfThisRow.second);
-                    }
-                }
-            }
-            //add self loops on the additional states (i.e., target and sink)
-            dtmcMatrixBuilder.addNextValue(targetState, targetState, storm::utility::one<ConstantType>());
-            dtmcMatrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one<ConstantType>());
-            mdpMatrixBuilder.newRowGroup(currentMdpRow);
-            mdpMatrixBuilder.addNextValue(currentMdpRow, targetState, storm::utility::one<ConstantType>());
-            ++currentMdpRow;
-            mdpMatrixBuilder.newRowGroup(currentMdpRow);
-            mdpMatrixBuilder.addNextValue(currentMdpRow, sinkState, storm::utility::one<ConstantType>());
-            ++currentMdpRow;
-            //The DTMC labeling
-            storm::models::sparse::StateLabeling dtmcStateLabeling(numStates);
-            storm::storage::BitVector initLabel(numStates, false);
-            initLabel.set(newStateIndexMap[this->initialState], true);
-            dtmcStateLabeling.addLabel("init", std::move(initLabel));
-            storm::storage::BitVector targetLabel(numStates, false);
-            targetLabel.set(numStates-2, true);
-            dtmcStateLabeling.addLabel("target", std::move(targetLabel));
-            storm::storage::BitVector sinkLabel(numStates, false);
-            sinkLabel.set(numStates-1, true);
-            dtmcStateLabeling.addLabel("sink", std::move(sinkLabel));
-            // The MDP labeling (is the same)
-            storm::models::sparse::StateLabeling mdpStateLabeling(dtmcStateLabeling);
-            // other ingredients
-            boost::optional<std::vector<ConstantType>> dtmcNoStateRewards, mdpNoStateRewards;
-            boost::optional<storm::storage::SparseMatrix<ConstantType>> dtmcNoTransitionRewards, mdpNoTransitionRewards;  
-            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> dtmcNoChoiceLabeling, mdpNoChoiceLabeling;            
-            // the final dtmc and mdp
-            sampleDtmc=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(dtmcMatrixBuilder.build(), std::move(dtmcStateLabeling), std::move(dtmcNoStateRewards), std::move(dtmcNoTransitionRewards), std::move(dtmcNoChoiceLabeling));
-            approxMdp=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(mdpMatrixBuilder.build(), std::move(mdpStateLabeling), std::move(mdpNoStateRewards), std::move(mdpNoTransitionRewards), std::move(mdpNoChoiceLabeling));
-            //build mapping vectors and the substitution vector. 
-            sampleDtmcMapping=std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>>();
-            sampleDtmcMapping.reserve(sampleDtmc->getTransitionMatrix().getEntryCount());
-            approxMdpMapping=std::vector<std::tuple<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&, std::size_t>>();
-            approxMdpMapping.reserve(approxMdp->getTransitionMatrix().getEntryCount());
-            approxMdpSubstitutions = std::vector<std::map<VariableType, TypeOfBound> >();
-            
-            currentMdpRow=0;
-            //go through rows:
-            for(storm::storage::sparse::state_type oldStateIndex : subsys){
-                //Get the index of the substitution for the current mdp row. (add it to approxMdpSubstitutions if we see this substitution the first time)
-                std::size_t substitutionIndex = std::find(approxMdpSubstitutions.begin(),approxMdpSubstitutions.end(), mdpRowSubstitutions[currentMdpRow]) - approxMdpSubstitutions.begin();
-                if(substitutionIndex==approxMdpSubstitutions.size()){
-                    approxMdpSubstitutions.push_back(mdpRowSubstitutions[currentMdpRow]);
-                }
-                typename storm::storage::SparseMatrix<ConstantType>::iterator dtmcEntryIt=sampleDtmc->getTransitionMatrix().getRow(newStateIndexMap[oldStateIndex]).begin();
-                typename storm::storage::SparseMatrix<ConstantType>::iterator mdpEntryIt=approxMdp->getTransitionMatrix().getRow(currentMdpRow).begin();
-                //go through columns to fill the dtmc and the first mdp row (of this group):
-                for(auto const& entry: transitions.getRow(oldStateIndex)){
-                    STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the DTMC matrix and the pDtmc Transitionmatrix do not match");
-                    STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the MDP matrix and the pDtmc Transitionmatrix do not match (" << newStateIndexMap[entry.getColumn()] << " vs " << mdpEntryIt->getColumn() << ")");
-                    sampleDtmcMapping.emplace_back(entry.getValue(),*dtmcEntryIt);
-                    approxMdpMapping.emplace_back(entry.getValue(),*mdpEntryIt, substitutionIndex);
-                    ++dtmcEntryIt;
-                    ++mdpEntryIt;
-                }
-                //transition to target state
-                if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ 
-                    STORM_LOG_THROW(targetState==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the DTMC matrix should match the target state but it does not.");
-                    STORM_LOG_THROW(targetState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the target state but it does not.");
-                    sampleDtmcMapping.emplace_back(oneStepProbs[oldStateIndex],*dtmcEntryIt);
-                    approxMdpMapping.emplace_back(oneStepProbs[oldStateIndex], *mdpEntryIt, substitutionIndex);
-                    ++dtmcEntryIt;
-                    ++mdpEntryIt;
-                }
-                //transition to sink state
-                if(!this->parametricTypeComparator.isZero(oneStepToSink[oldStateIndex])){ 
-                    STORM_LOG_THROW(sinkState==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the DTMC matrix should match the sink state but it does not.");
-                    STORM_LOG_THROW(sinkState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the sink state but it does not.");
-                    sampleDtmcMapping.emplace_back(oneStepToSink[oldStateIndex],*dtmcEntryIt);
-                    approxMdpMapping.emplace_back(oneStepToSink[oldStateIndex], *mdpEntryIt, substitutionIndex);
-                    ++dtmcEntryIt;
-                    ++mdpEntryIt;
-                }
-                //fill in the remaining rows of the mdp
-                storm::storage::sparse::state_type startOfNextRowGroup = currentMdpRow + approxMdp->getTransitionMatrix().getRowGroupSize(newStateIndexMap[oldStateIndex]);
-                for(++currentMdpRow;currentMdpRow<startOfNextRowGroup; ++currentMdpRow){
-                    //Get the new substitution index
-                    substitutionIndex = std::find(approxMdpSubstitutions.begin(),approxMdpSubstitutions.end(), mdpRowSubstitutions[currentMdpRow]) - approxMdpSubstitutions.begin();
-                    if(substitutionIndex==approxMdpSubstitutions.size()){
-                        approxMdpSubstitutions.push_back(mdpRowSubstitutions[currentMdpRow]);
-                    }
-                    //go through columns
-                    mdpEntryIt=approxMdp->getTransitionMatrix().getRow(currentMdpRow).begin();
-                    for(auto const& entry: transitions.getRow(oldStateIndex)){
-                        STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the MDP matrix and the pDtmc Transitionmatrix do not match");
-                        approxMdpMapping.emplace_back(entry.getValue(),*mdpEntryIt, substitutionIndex);
-                        ++mdpEntryIt;
-                    }
-                    //transition to target state
-                    if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ 
-                        STORM_LOG_THROW(targetState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the target state but it does not.");
-                        approxMdpMapping.emplace_back(oneStepProbs[oldStateIndex], *mdpEntryIt, substitutionIndex);
-                        ++mdpEntryIt;
-                    }
-                    //transition to sink state
-                    if(!this->parametricTypeComparator.isZero(oneStepToSink[oldStateIndex])){ 
-                        STORM_LOG_THROW(sinkState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the sink state but it does not.");
-                        approxMdpMapping.emplace_back(oneStepToSink[oldStateIndex], *mdpEntryIt, substitutionIndex);
-                        ++mdpEntryIt;
-                    }
-                }
-            }
-        }
-        */
-        template<typename ParametricType, typename ConstantType>
-        ParametricType SparseDtmcRegionModelChecker<ParametricType, ConstantType>::getReachProbFunction() {
-            if(!this->isReachProbFunctionComputed){
-                std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionStart = std::chrono::high_resolution_clock::now();
-                //get the one step probabilities and the transition matrix of the simplified model without target/sink state
-                //TODO check if this takes long... we could also store the oneSteps while creating the simplified model. Or(?) we keep the matrix the way it is and give the target state one step probability 1.
-                storm::storage::SparseMatrix<ParametricType> backwardTransitions(this->simplifiedModel->getBackwardTransitions());
-                std::vector<ParametricType> oneStepProbabilities(this->simplifiedModel->getNumberOfStates()-2, storm::utility::zero<ParametricType>());
-                for(auto const& entry : backwardTransitions.getRow(*(this->simplifiedModel->getStates("target").begin()))){
-                    if(entry.getColumn()<oneStepProbabilities.size()){
-                        oneStepProbabilities[entry.getColumn()]=entry.getValue();
-                    } //else case: only holds for the entry that corresponds to the selfloop on the target state..
-                }
-                storm::storage::BitVector maybeStates=~(this->simplifiedModel->getStates("target") | this->simplifiedModel->getStates("sink"));
-                backwardTransitions=backwardTransitions.getSubmatrix(false,maybeStates,maybeStates);
-                storm::storage::SparseMatrix<ParametricType> forwardTransitions=this->simplifiedModel->getTransitionMatrix().getSubmatrix(false,maybeStates,maybeStates);
-                //now compute the functions using methods from elimination model checker
-                storm::storage::BitVector newInitialStates = this->simplifiedModel->getInitialStates() % maybeStates;
-                storm::storage::BitVector phiStates(this->simplifiedModel->getNumberOfStates(), true);
-                boost::optional<std::vector<ParametricType>> noStateRewards;
-                std::vector<std::size_t> statePriorities = this->eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities);
-                this->reachProbFunction=this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, this->simplifiedModel->getStates("target"), noStateRewards, statePriorities);
-                std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionEnd = std::chrono::high_resolution_clock::now();
-                std::string funcStr = " (/ " +
-                                    this->reachProbFunction.nominator().toString(false, true) + " " +
-                                    this->reachProbFunction.denominator().toString(false, true) +
-                                " )";
-                std::cout << std::endl <<"the resulting reach prob function is " << std::endl << funcStr << std::endl << std::endl;
-                STORM_LOG_DEBUG("Computed reachProbFunction");
-                this->isReachProbFunctionComputed=true;
-                this->timeComputeReachProbFunction = timeComputeReachProbFunctionEnd-timeComputeReachProbFunctionStart;
-            }
-            return this->reachProbFunction;
-        }
-
-        
-
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, const ParametricType& reachProbFunc, const storm::logic::ProbabilityOperatorFormula& formula) {
-                    
-            storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions..
-            solver = std::shared_ptr<storm::solver::Smt2SmtSolver>(new storm::solver::Smt2SmtSolver(manager, true));
-            
-            ParametricType bound= storm::utility::regions::convertNumber<double, ParametricType>(this->probabilityOperatorFormula->getBound());
-            
-            // To prove that the property is satisfied in the initial state for all parameters,
-            // we ask the solver whether the negation of the property is satisfiable and invert the answer.
-            // In this case, assert that this variable is true:
-            VariableType proveAllSatVar=storm::utility::regions::getNewVariable<VariableType>("storm_proveAllSat", storm::utility::regions::VariableSort::VS_BOOL);
-            
-            //Example:
-            //Property:    P<=p [ F 'target' ] holds iff...
-            // f(x)         <= p
-            // Hence: If  f(x) > p is unsat, the property is satisfied for all parameters.
-            
-            storm::logic::ComparisonType proveAllSatRel; //the relation from the property needs to be inverted
-            switch (this->probabilityOperatorFormula->getComparisonType()) {
-                case storm::logic::ComparisonType::Greater:
-                    proveAllSatRel=storm::logic::ComparisonType::LessEqual;
-                    break;
-                case storm::logic::ComparisonType::GreaterEqual:
-                    proveAllSatRel=storm::logic::ComparisonType::Less;
-                    break;
-                case storm::logic::ComparisonType::Less:
-                    proveAllSatRel=storm::logic::ComparisonType::GreaterEqual;
-                    break;
-                case storm::logic::ComparisonType::LessEqual:
-                    proveAllSatRel=storm::logic::ComparisonType::Greater;
-                    break;
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-            }
-            storm::utility::regions::addGuardedConstraintToSmtSolver(solver, proveAllSatVar, reachProbFunc, proveAllSatRel, bound);
-            
-            // To prove that the property is violated in the initial state for all parameters,
-            // we ask the solver whether the the property is satisfiable and invert the answer.
-            // In this case, assert that this variable is true:
-            VariableType proveAllViolatedVar=storm::utility::regions::getNewVariable<VariableType>("storm_proveAllViolated", storm::utility::regions::VariableSort::VS_BOOL);        
-            
-            //Example:
-            //Property:    P<=p [ F 'target' ] holds iff...
-            // f(x)         <= p
-            // Hence: If f(x)  <= p is unsat, the property is violated for all parameters. 
-            storm::logic::ComparisonType proveAllViolatedRel = this->probabilityOperatorFormula->getComparisonType();
-            storm::utility::regions::addGuardedConstraintToSmtSolver(solver, proveAllViolatedVar, reachProbFunc, proveAllViolatedRel, bound);          
-        }
-
-
-
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkRegion(ParameterRegion& region) {
-            std::chrono::high_resolution_clock::time_point timeCheckRegionStart = std::chrono::high_resolution_clock::now();
-            ++this->numOfCheckedRegions;
-            
-            STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to analyze a region although no property has been specified" );
-            STORM_LOG_DEBUG("Analyzing the region " << region.toString());
-           // std::cout << "Analyzing the region " << region.toString() << std::endl;
-            
-            //switches for the different steps.
-            bool doApproximation=this->hasOnlyLinearFunctions; // this approach is only correct if the model has only linear functions
-            bool doSampling=true;
-            bool doSubsystemSmt=false; //this->hasOnlyLinearFunctions;  // this approach is only correct if the model has only linear functions
-            bool doFullSmt=false;
-            
-            std::chrono::high_resolution_clock::time_point timeApproximationStart = std::chrono::high_resolution_clock::now();
-            std::vector<ConstantType> lowerBounds;
-            std::vector<ConstantType> upperBounds;
-            if(doApproximation){
-                STORM_LOG_DEBUG("Checking approximative probabilities...");
-          //      std::cout << "  Checking approximative probabilities..." << std::endl;
-                if(checkApproximativeProbabilities(region, lowerBounds, upperBounds)){
-                    ++this->numOfRegionsSolvedThroughApproximation;
-                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through approximation.");
-                    doSampling=false;
-                    doSubsystemSmt=false;
-                    doFullSmt=false;
-                }
-            }
-            std::chrono::high_resolution_clock::time_point timeApproximationEnd = std::chrono::high_resolution_clock::now();
-            
-            std::chrono::high_resolution_clock::time_point timeSamplingStart = std::chrono::high_resolution_clock::now();
-            if(doSampling){
-                STORM_LOG_DEBUG("Checking sample points...");
-        //        std::cout << "  Checking sample points..." << std::endl;
-                if(checkSamplePoints(region)){
-                    ++this->numOfRegionsSolvedThroughSampling;
-                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through sampling.");
-                    doApproximation=false;
-                    doSubsystemSmt=false;
-                    doFullSmt=false;
-                }
-            }
-            std::chrono::high_resolution_clock::time_point timeSamplingEnd = std::chrono::high_resolution_clock::now();
-            
-            std::chrono::high_resolution_clock::time_point timeSubsystemSmtStart = std::chrono::high_resolution_clock::now();
-            if(doSubsystemSmt){
-                STORM_LOG_WARN("SubsystemSmt approach not yet implemented");
-            }
-            std::chrono::high_resolution_clock::time_point timeSubsystemSmtEnd = std::chrono::high_resolution_clock::now();
-            
-            std::chrono::high_resolution_clock::time_point timeFullSmtStart = std::chrono::high_resolution_clock::now();
-            if(doFullSmt){
-                STORM_LOG_DEBUG("Checking with Smt Solving...");
-     //           std::cout << "  Checking with Smt Solving..." << std::endl;
-                if(checkFullSmt(region)){
-                    ++this->numOfRegionsSolvedThroughFullSmt;
-                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through Smt Solving.");
-                    doApproximation=false;
-                    doSampling=false;
-                    doSubsystemSmt=false;
-                }
-            }
-            std::chrono::high_resolution_clock::time_point timeFullSmtEnd = std::chrono::high_resolution_clock::now();
-            
-            
-            //some information for statistics...
-            std::chrono::high_resolution_clock::time_point timeCheckRegionEnd = std::chrono::high_resolution_clock::now();
-            this->timeCheckRegion += timeCheckRegionEnd-timeCheckRegionStart;
-            this->timeSampling += timeSamplingEnd - timeSamplingStart;
-            this->timeApproximation += timeApproximationEnd - timeApproximationStart;
-            this->timeSubsystemSmt += timeSubsystemSmtEnd - timeSubsystemSmtStart;
-            this->timeFullSmt += timeFullSmtEnd - timeFullSmtStart;
-            switch(region.getCheckResult()){
-                case RegionCheckResult::EXISTSBOTH:
-                    ++this->numOfRegionsExistsBoth;
-                    break;
-                case RegionCheckResult::ALLSAT:
-                    ++this->numOfRegionsAllSat;
-                    break;
-                case RegionCheckResult::ALLVIOLATED:
-                    ++this->numOfRegionsAllViolated;
-                    break;
-                default:
-                    break;
-            }
-            
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkRegions(std::vector<ParameterRegion>& regions) {
-            STORM_LOG_DEBUG("Checking " << regions.size() << "regions.");
-            std::cout << "Checking " << regions.size() << "regions. Progress: ";
-            std::cout.flush();
-                    
-            uint_fast64_t progress=0;
-            uint_fast64_t checkedRegions=0;
-            for(auto& region : regions){
-                this->checkRegion(region);
-                if((checkedRegions++)*10/regions.size()==progress){
-                    std::cout << progress++;
-                    std::cout.flush();
-                }
-            }
-            std::cout << std::endl;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkSamplePoints(ParameterRegion& region) {
-            
-            auto samplingPoints = region.getVerticesOfRegion(region.getVariables()); //test the 4 corner points
-            for (auto const& point : samplingPoints){
-                if(checkPoint(region, point, false)){
-                    return true;
-                }            
-            }
-            return false;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkPoint(ParameterRegion& region, std::map<VariableType, CoefficientType>const& point, bool viaReachProbFunction) {
-            bool valueInBoundOfFormula;
-            if(viaReachProbFunction){
-                valueInBoundOfFormula = this->valueIsInBoundOfFormula(storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(getReachProbFunction(), point));
-            }
-            else{
-                this->samplingModel->instantiate(point);
-                valueInBoundOfFormula=this->valueIsInBoundOfFormula(this->samplingModel->computeReachabilityProbabilities()[*this->samplingModel->getModel()->getInitialStates().begin()]);
-            }
-                
-            if(valueInBoundOfFormula){
-                if (region.getCheckResult()!=RegionCheckResult::EXISTSSAT){
-                    region.setSatPoint(point);
-                    if(region.getCheckResult()==RegionCheckResult::EXISTSVIOLATED){
-                        region.setCheckResult(RegionCheckResult::EXISTSBOTH);
-                        return true;
-                    }
-                    region.setCheckResult(RegionCheckResult::EXISTSSAT);
-                }
-            }
-            else{
-                if (region.getCheckResult()!=RegionCheckResult::EXISTSVIOLATED){
-                    region.setViolatedPoint(point);
-                    if(region.getCheckResult()==RegionCheckResult::EXISTSSAT){
-                        region.setCheckResult(RegionCheckResult::EXISTSBOTH);
-                        return true;
-                    }
-                    region.setCheckResult(RegionCheckResult::EXISTSVIOLATED);
-                }
-            }
-            return false;
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkApproximativeProbabilities(ParameterRegion& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds) {
-                    
-            STORM_LOG_THROW(this->hasOnlyLinearFunctions, storm::exceptions::UnexpectedException, "Tried to perform approximative method (only applicable if all functions are linear), but there are nonlinear functions.");
-            std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now();
-            this->approximationModel->instantiate(region);
-            std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now();
-            this->timeMDPBuild += timeMDPBuildEnd-timeMDPBuildStart;
-            
-            if (region.getCheckResult()==RegionCheckResult::UNKNOWN){
-                //Sampling to get a hint of whether we should try to prove ALLSAT or ALLVIOLATED
-                checkPoint(region,region.getLowerBounds(), false);
-            }
-            
-            // Decide whether we should compute lower or upper bounds first.
-            storm::logic::OptimalityType firstOpType;
-            switch (this->probabilityOperatorFormula->getComparisonType()) {
-                case storm::logic::ComparisonType::Greater:
-                case storm::logic::ComparisonType::GreaterEqual:
-                    switch (region.getCheckResult()){
-                        case RegionCheckResult::EXISTSSAT:
-                        case RegionCheckResult::UNKNOWN: 
-                            firstOpType=storm::logic::OptimalityType::Minimize;
-                            break;
-                        case RegionCheckResult::EXISTSVIOLATED:
-                            firstOpType=storm::logic::OptimalityType::Maximize;
-                            break;
-                        default:
-                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The checkresult of the current region should not be conclusive, i.e. it should be either EXISTSSAT or EXISTSVIOLATED or UNKNOWN");
-                    }
-                    break;
-                case storm::logic::ComparisonType::Less:
-                case storm::logic::ComparisonType::LessEqual:
-                    switch (region.getCheckResult()){
-                        case RegionCheckResult::EXISTSSAT:
-                        case RegionCheckResult::UNKNOWN: 
-                            firstOpType=storm::logic::OptimalityType::Maximize;
-                            break;
-                        case RegionCheckResult::EXISTSVIOLATED:
-                            firstOpType=storm::logic::OptimalityType::Minimize;
-                            break;
-                        default:
-                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The checkresult of the current region should not be conclusive, i.e. it should be either EXISTSSAT or EXISTSVIOLATED or UNKNOWN");
-                    }
-                    break;
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-            }
-            
-            //one iteration for each opType \in {Maximize, Minimize}
-            //However, the second iteration might not be performed if the first one proved ALLSAT or ALLVIOLATED
-            storm::logic::OptimalityType opType=firstOpType;
-            for(int i=1; i<=2; ++i){   
-                
-                //check if the approximation yielded a conclusive result
-                if(opType==storm::logic::OptimalityType::Maximize){
-                    upperBounds = this->approximationModel->computeReachabilityProbabilities(opType);
-                    if(valueIsInBoundOfFormula(upperBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
-                        if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Less) || 
-                                (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::LessEqual)){
-                            region.setCheckResult(RegionCheckResult::ALLSAT);
-                            return true;
-                        }
-                    }
-                    else{
-                        if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Greater) || 
-                                (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::GreaterEqual)){
-                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
-                            return true;
-                        }
-                    }
-                    //flip the optype for the next iteration
-                    opType=storm::logic::OptimalityType::Minimize;
-                }
-                else if(opType==storm::logic::OptimalityType::Minimize){
-                    lowerBounds = this->approximationModel->computeReachabilityProbabilities(opType);
-                    if(valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
-                        if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Greater) || 
-                                (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::GreaterEqual)){
-                            region.setCheckResult(RegionCheckResult::ALLSAT);
-                            return true;
-                        }
-                    }
-                    else{
-                        if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Less) || 
-                                (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::LessEqual)){
-                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
-                            return true;
-                        }
-                    }                
-                    //flip the optype for the next iteration
-                    opType=storm::logic::OptimalityType::Maximize;
-                }
-            }
-            //if we reach this point than the result is still inconclusive.
-            return false;            
-        }
-        
-        /*
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::buildMdpForApproximation(const ParameterRegion& region) {
-            
-            //instantiate the substitutions for the given region
-            std::vector<std::map<VariableType, CoefficientType>> substitutions(this->approxMdpSubstitutions.size());
-            for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->approxMdpSubstitutions.size(); ++substitutionIndex){
-                for(std::pair<VariableType, TypeOfBound> const& sub : this->approxMdpSubstitutions[substitutionIndex]){
-                    switch(sub.second){
-                        case TypeOfBound::LOWER:
-                            substitutions[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first)));
-                            break;
-                        case TypeOfBound::UPPER:
-                            substitutions[substitutionIndex].insert(std::make_pair(sub.first, region.getUpperBound(sub.first)));
-                            break;
-                        default:
-                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound");
-                    }
-                }
-            }
-            
-            std::vector<std::pair<typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, ConstantType>> entryVector;
-            entryVector.reserve(this->approxMdpMapping.size());
-            //now put the values into the mdp matrix
-            for( std::tuple<ParametricType, typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, size_t>& mappingTriple : this->approxMdpMapping){
-                entryVector.emplace_back(std::get<1>(mappingTriple), storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
-                                                storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)])));
-            }
-            
-            
-            
-            std::chrono::high_resolution_clock::time_point timeMDPInsertingStart = std::chrono::high_resolution_clock::now();
-            for(std::pair<typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, ConstantType>& entry : entryVector ){
-                entry.first.setValue(entry.second);
-            }
-            std::chrono::high_resolution_clock::time_point timeMDPInsertingEnd = std::chrono::high_resolution_clock::now();
-            this->timeMDPInserting += timeMDPInsertingEnd-timeMDPInsertingStart;
-    
-        }
-*/        
-      
-
-        
-        
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkFullSmt(ParameterRegion& region) {
-            if (region.getCheckResult()==RegionCheckResult::UNKNOWN){
-                //Sampling needs to be done (on a single point)
-                checkPoint(region,region.getLowerBounds(), true);
-            }
-            
-            if(this->smtSolver==nullptr){
-                initializeSMTSolver(this->smtSolver, getReachProbFunction(),*this->probabilityOperatorFormula);
-            }
-            
-            this->smtSolver->push();
-            
-            //add constraints for the region
-            for(auto const& variable : region.getVariables()) {
-                storm::utility::regions::addParameterBoundsToSmtSolver(this->smtSolver, variable, storm::logic::ComparisonType::GreaterEqual, region.getLowerBound(variable));
-                storm::utility::regions::addParameterBoundsToSmtSolver(this->smtSolver, variable, storm::logic::ComparisonType::LessEqual, region.getUpperBound(variable));
-            }
-            
-            //add constraint that states what we want to prove            
-            VariableType proveAllSatVar=storm::utility::regions::getVariableFromString<VariableType>("storm_proveAllSat");     
-            VariableType proveAllViolatedVar=storm::utility::regions::getVariableFromString<VariableType>("storm_proveAllViolated");     
-            switch(region.getCheckResult()){
-                case RegionCheckResult::EXISTSBOTH:
-                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::EXISTSBOTH), "checkFullSmt invoked although the result is already clear (EXISTSBOTH). Will validate this now...");
-                case RegionCheckResult::ALLSAT:
-                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::ALLSAT), "checkFullSmt invoked although the result is already clear (ALLSAT). Will validate this now...");
-                case RegionCheckResult::EXISTSSAT:
-                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllSatVar, true);
-                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllViolatedVar, false);
-                    break;
-                case RegionCheckResult::ALLVIOLATED:
-                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::ALLVIOLATED), "checkFullSmt invoked although the result is already clear (ALLVIOLATED). Will validate this now...");
-                case RegionCheckResult::EXISTSVIOLATED:
-                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllSatVar, false);
-                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllViolatedVar, true);
-                    break;
-                default:
-                STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not handle the current region CheckResult: " << region.checkResultToString());
-            }
-            
-            storm::solver::SmtSolver::CheckResult solverResult= this->smtSolver->check();
-            this->smtSolver->pop();
-            
-            switch(solverResult){
-                case storm::solver::SmtSolver::CheckResult::Sat:
-                    switch(region.getCheckResult()){
-                        case RegionCheckResult::EXISTSSAT:
-                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
-                            //There is also a violated point
-                            STORM_LOG_WARN("Extracting a violated point from the smt solver is not yet implemented!");
-                            break;
-                        case RegionCheckResult::EXISTSVIOLATED:
-                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
-                            //There is also a sat point
-                            STORM_LOG_WARN("Extracting a sat point from the smt solver is not yet implemented!");
-                            break;
-                        case RegionCheckResult::EXISTSBOTH:
-                            //That was expected
-                            STORM_LOG_WARN("result EXISTSBOTH Validated!");
-                            break;
-                        default:
-                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The solver gave an unexpected result (sat)");
-                    }
-                    return true;
-                case storm::solver::SmtSolver::CheckResult::Unsat:
-                    switch(region.getCheckResult()){
-                        case RegionCheckResult::EXISTSSAT:
-                            region.setCheckResult(RegionCheckResult::ALLSAT);
-                            break;
-                        case RegionCheckResult::EXISTSVIOLATED:
-                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
-                            break;
-                        case RegionCheckResult::ALLSAT:
-                            //That was expected...
-                            STORM_LOG_WARN("result ALLSAT Validated!");
-                            break;
-                        case RegionCheckResult::ALLVIOLATED:
-                            //That was expected...
-                            STORM_LOG_WARN("result ALLVIOLATED Validated!");
-                            break;
-                        default:
-                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The solver gave an unexpected result (unsat)");
-                    }
-                    return true;
-                case storm::solver::SmtSolver::CheckResult::Unknown:
-                default:
-                    STORM_LOG_WARN("The SMT solver was not able to compute a result for this region. (Timeout? Memout?)");
-                    if(this->smtSolver->isNeedsRestart()){
-                        initializeSMTSolver(this->smtSolver,getReachProbFunction(), *this->probabilityOperatorFormula);
-                    }
-                    return false;
-            }
-        }
-
-             template<typename ParametricType, typename ConstantType>
-        template<typename ValueType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::valueIsInBoundOfFormula(ValueType value){
-            STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to compare a value to the bound of a formula, but no formula specified.");
-            double valueAsDouble = storm::utility::regions::convertNumber<ValueType, double>(value);
-        //    std::cout << "checked whether value: " << value << " (= " << valueAsDouble << " double ) is in bound of formula: " << *this->probabilityOperatorFormula << std::endl;
-        //    std::cout << "     (valueAsDouble >= this->probabilityOperatorFormula->getBound()) is " << (valueAsDouble >= this->probabilityOperatorFormula->getBound()) << std::endl;
-            switch (this->probabilityOperatorFormula->getComparisonType()) {
-                case storm::logic::ComparisonType::Greater:
-                    return (valueAsDouble > this->probabilityOperatorFormula->getBound());
-                case storm::logic::ComparisonType::GreaterEqual:
-                    return (valueAsDouble >= this->probabilityOperatorFormula->getBound());
-                case storm::logic::ComparisonType::Less:
-                    return (valueAsDouble < this->probabilityOperatorFormula->getBound());
-                case storm::logic::ComparisonType::LessEqual:
-                    return (valueAsDouble <= this->probabilityOperatorFormula->getBound());
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-            }
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::printStatisticsToStream(std::ostream& outstream) {
-            
-            if(this->probabilityOperatorFormula==nullptr){
-                outstream << "Statistic Region Model Checker Statistics Error: No formula specified." << std::endl; 
-                return;
-            }
-            
-            std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timePreprocessing);
-            std::chrono::milliseconds timeInitialStateEliminationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeInitialStateElimination);
-            std::chrono::milliseconds timeComputeReachProbFunctionInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeComputeReachProbFunction);
-            std::chrono::milliseconds timeCheckRegionInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeCheckRegion);
-            std::chrono::milliseconds timeSammplingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeSampling);
-            std::chrono::milliseconds timeApproximationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeApproximation);
-            std::chrono::milliseconds timeMDPBuildInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeMDPBuild);
-            std::chrono::milliseconds timeFullSmtInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeFullSmt);
-            
-            std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ...
-            std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeOverall);
-            
-            uint_fast64_t numOfSolvedRegions= this->numOfRegionsExistsBoth + this->numOfRegionsAllSat + this->numOfRegionsAllViolated;
-            
-            outstream << std::endl << "Statistics Region Model Checker Statistics:" << std::endl;
-            outstream << "-----------------------------------------------" << std::endl;
-            outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl;
-            outstream << "Simplified model: " << this->simplifiedModel->getNumberOfStates() << " states, " << this->simplifiedModel->getNumberOfTransitions() << " transitions" << std::endl;
-            outstream << "Formula: " << *this->probabilityOperatorFormula << std::endl;
-            outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl;
-            outstream << "Number of checked regions: " << this->numOfCheckedRegions << std::endl;
-            outstream << "  Number of solved regions:  " <<  numOfSolvedRegions << "(" << numOfSolvedRegions*100/this->numOfCheckedRegions << "%)" <<  std::endl;
-            outstream << "    AllSat:      " <<  this->numOfRegionsAllSat << "(" << this->numOfRegionsAllSat*100/this->numOfCheckedRegions << "%)" <<  std::endl;
-            outstream << "    AllViolated: " <<  this->numOfRegionsAllViolated << "(" << this->numOfRegionsAllViolated*100/this->numOfCheckedRegions << "%)" <<  std::endl;
-            outstream << "    ExistsBoth:  " <<  this->numOfRegionsExistsBoth << "(" << this->numOfRegionsExistsBoth*100/this->numOfCheckedRegions << "%)" <<  std::endl;
-            outstream << "    Unsolved:    " <<  this->numOfCheckedRegions - numOfSolvedRegions << "(" << (this->numOfCheckedRegions - numOfSolvedRegions)*100/this->numOfCheckedRegions << "%)" <<  std::endl;
-            outstream << "  --  Note: %-numbers are relative to the NUMBER of regions, not the size of their area --" <<  std::endl;
-            outstream << "  " << this->numOfRegionsSolvedThroughSampling << " regions solved through Sampling" << std::endl;
-            outstream << "  " << this->numOfRegionsSolvedThroughApproximation << " regions solved through Approximation" << std::endl;
-            outstream << "  " << this->numOfRegionsSolvedThroughSubsystemSmt << " regions solved through SubsystemSmt" << std::endl;
-            outstream << "  " << this->numOfRegionsSolvedThroughFullSmt << " regions solved through FullSmt" << std::endl;
-            outstream << std::endl;
-            outstream << "Running times:" << std::endl;
-            outstream << "  " << timeOverallInMilliseconds.count() << "ms overall" << std::endl;
-            outstream << "  " << timePreprocessingInMilliseconds.count() << "ms Preprocessing including... " << std::endl;
-            outstream << "    " << timeInitialStateEliminationInMilliseconds.count() << "ms Initial state elimination of const transitions" << std::endl;
-            outstream << "  " << timeComputeReachProbFunctionInMilliseconds.count() << "ms to compute the reachability probability function" << std::endl;
-            outstream << "  " << timeCheckRegionInMilliseconds.count() << "ms Region Check including... " << std::endl;
-            outstream << "    " << timeSammplingInMilliseconds.count() << "ms Sampling " << std::endl;
-            outstream << "    " << timeApproximationInMilliseconds.count() << "ms Approximation including... " << std::endl;
-            outstream << "      " << timeMDPBuildInMilliseconds.count() << "ms to build the MDP" << std::endl;
-            outstream << "    " << timeFullSmtInMilliseconds.count() << "ms Full Smt solving" << std::endl;
-            outstream << "-----------------------------------------------" << std::endl;
-            
-        }
-   
-
-
-
-#ifdef STORM_HAVE_CARL
-        
-        //DELETEME
-        template<>
-        std::pair<storm::storage::SparseMatrix<double>,std::vector<boost::container::flat_set<uint_fast64_t>>> SparseDtmcRegionModelChecker<storm::RationalFunction, double>::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector<std::map<storm::Variable, storm::RationalFunction::CoeffType>> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector<storm::RationalFunction> const& oneStepProbabilities, bool addSelfLoops) const {
-            
-            //Check if the arguments are as expected
-            STORM_LOG_THROW((std::is_same<storm::RationalFunction::CoeffType, cln::cl_RA>::value), storm::exceptions::IllegalArgumentException, "Unexpected Type of Coefficients");
-            STORM_LOG_THROW(filter.size()==matrix.getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the filter");
-            STORM_LOG_THROW(oneStepProbabilities.empty() || oneStepProbabilities.size()==matrix.getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the oneStepProbabilities");
-            
-            //get a mapping from old state indices to the new ones
-            std::vector<storm::storage::sparse::state_type> newStateIndexMap(matrix.getNumberOfRows(), matrix.getNumberOfRows()); //initialize with some illegal index to easily check if a transition leads to an unselected state
-            storm::storage::sparse::state_type newStateIndex=0;
-            for(auto const& state : filter){
-                newStateIndexMap[state]=newStateIndex;
-                ++newStateIndex;
-            }
-            storm::storage::sparse::state_type numStates=filter.getNumberOfSetBits();
-            STORM_LOG_ASSERT(newStateIndex==numStates, "unexpected number of new states");
-            storm::storage::sparse::state_type targetState =0;
-            storm::storage::sparse::state_type sinkState=0;
-            if(!oneStepProbabilities.empty()){
-                targetState=numStates;
-                ++numStates;
-            }
-            if(addSinkState){
-                sinkState=numStates;
-                ++numStates;
-            }
-            //todo rows (i.e. the first parameter) should be numStates*substitutions.size ?
-            storm::storage::SparseMatrixBuilder<double> matrixBuilder(0, numStates, 0, true, true, numStates);
-            std::vector<boost::container::flat_set<uint_fast64_t>> choiceLabeling;
-            //fill in the data row by row            
-            storm::storage::sparse::state_type currentRow=0;
-            for(auto const& oldStateIndex : filter){
-                matrixBuilder.newRowGroup(currentRow);
-                for(size_t substitutionIndex=0; substitutionIndex< substitutions.size(); ++substitutionIndex){
-                    double missingProbability = 1.0;
-                    if(matrix.getRow(oldStateIndex).empty()){ //just add the selfloop if there is no transition
-                        if(addSelfLoops){
-                            matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero<double>());
-                        }
-                    }
-                    else{
-                        FlexibleMatrix::const_iterator entry = matrix.getRow(oldStateIndex).begin();
-                        for(; entry<matrix.getRow(oldStateIndex).end() && entry->getColumn()<oldStateIndex; ++entry){ //insert until we come to the diagonal entry
-                            double value = cln::double_approx(entry->getValue().evaluate(substitutions[substitutionIndex]));
-                            missingProbability-=value;
-                            storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()];
-                            STORM_LOG_THROW(column<numStates, storm::exceptions::IllegalArgumentException, "Illegal filter: Selected a state that has a transition to an unselected state.");
-                            matrixBuilder.addNextValue(currentRow, column, value);
-                        }
-                        if(addSelfLoops && entry->getColumn()!=oldStateIndex){ //maybe add a zero valued diagonal entry
-                            matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero<double>());
-                        }
-                        for(; entry < matrix.getRow(oldStateIndex).end(); ++entry){ //insert the rest
-                            double value = cln::double_approx(entry->getValue().evaluate(substitutions[substitutionIndex]));
-                            missingProbability-=value;
-                            storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()];
-                            STORM_LOG_THROW(column<numStates, storm::exceptions::IllegalArgumentException, "Illegal filter: Selected a state that has a transition to an unselected state.");
-                            matrixBuilder.addNextValue(currentRow, column, value);
-                        }
-                    }
-                    if(!oneStepProbabilities.empty() && !oneStepProbabilities[oldStateIndex].isZero()){ //transition to target state
-                        double value = cln::double_approx(oneStepProbabilities[oldStateIndex].evaluate(substitutions[substitutionIndex]));
-                        missingProbability-=value;
-                        matrixBuilder.addNextValue(currentRow, targetState, value);
-                    }
-                    storm::utility::ConstantsComparator<double> doubleComperator;
-                    if(addSinkState && !doubleComperator.isZero(missingProbability)){ //transition to sink state
-                        STORM_LOG_ASSERT(missingProbability> -storm::settings::generalSettings().getPrecision(), "The missing probability is negative.");
-                        matrixBuilder.addNextValue(currentRow, sinkState, missingProbability);
-                    }
-                    boost::container::flat_set<uint_fast64_t> label;
-                    label.insert(substitutionIndex);
-                    choiceLabeling.emplace_back(label);
-                    ++currentRow;
-                }
-            }
-            //finally, add self loops on the additional states (i.e., target and sink)
-            boost::container::flat_set<uint_fast64_t> labelAll;
-            labelAll.insert(substitutions.size());
-            if (!oneStepProbabilities.empty()){
-                matrixBuilder.newRowGroup(currentRow);
-                matrixBuilder.addNextValue(currentRow, targetState, storm::utility::one<double>());
-                choiceLabeling.emplace_back(labelAll);
-                ++currentRow;
-            }
-            
-            if (addSinkState){
-                matrixBuilder.newRowGroup(currentRow);
-                matrixBuilder.addNextValue(currentRow, sinkState, storm::utility::one<double>());
-                choiceLabeling.emplace_back(labelAll);
-                ++currentRow;
-            }
-       
-            return std::pair<storm::storage::SparseMatrix<double>, std::vector<boost::container::flat_set<uint_fast64_t>>>(matrixBuilder.build(), std::move(choiceLabeling));
-        }
-        
-        //DELETEME
-        template<typename ParametricType, typename ConstantType>
-        std::pair<storm::storage::SparseMatrix<double>,std::vector<boost::container::flat_set<uint_fast64_t>>> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector<std::map<storm::Variable, storm::RationalFunction::CoeffType>> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector<ParametricType> const& oneStepProbabilities, bool addSelfLoops) const{
-            STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Instantiation of flexible matrix is not supported for this type");
-        }
-
-
-      //DELETEME
-        template<>
-        void SparseDtmcRegionModelChecker<storm::RationalFunction, double>::eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector<storm::RationalFunction>& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix<storm::RationalFunction> const& forwardTransitions, boost::optional<std::vector<std::size_t>> const& statePriorities){
-    
-            if(true){ // eliminate all states with constant outgoing transitions
-                storm::storage::BitVector statesToEliminate = ~initialstates;
-                //todo: ordering of states important?
-                boost::optional<std::vector<storm::RationalFunction>> missingStateRewards;
-                for (auto const& state : statesToEliminate) {
-                    bool onlyConstantOutgoingTransitions=true;
-                    for(auto const& entry : flexibleMatrix.getRow(state)){
-                        if(!entry.getValue().isConstant()){
-                            onlyConstantOutgoingTransitions=false;
-                            break;
-                        }
-                    }
-                    if(onlyConstantOutgoingTransitions){
-                        eliminationModelChecker.eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
-                        subsystem.set(state,false);
-                    }
-                }
-
-                //Note: we could also "eliminate" the initial state to get rid of its selfloop
-            }
-            else if(false){ //eliminate all states with standard state elimination
-                boost::optional<std::vector<storm::RationalFunction>> missingStateRewards;
-                
-                storm::storage::BitVector statesToEliminate = ~initialstates;
-                std::vector<storm::storage::sparse::state_type> states(statesToEliminate.begin(), statesToEliminate.end());
-                
-                if (statePriorities) {
-                    std::sort(states.begin(), states.end(), [&statePriorities] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statePriorities.get()[a] < statePriorities.get()[b]; });
-                }
-                
-                STORM_LOG_DEBUG("Eliminating " << states.size() << " states using the state elimination technique." << std::endl);
-                for (auto const& state : states) {
-                    eliminationModelChecker.eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
-                }
-                subsystem=~statesToEliminate;
-                
-            }
-            else if(false){ //hybrid method
-                boost::optional<std::vector<storm::RationalFunction>> missingStateRewards;
-                storm::storage::BitVector statesToEliminate = ~initialstates;
-                uint_fast64_t maximalDepth =0;
-                std::vector<storm::storage::sparse::state_type> entryStateQueue;
-                STORM_LOG_DEBUG("Eliminating " << statesToEliminate.size() << " states using the hybrid elimination technique." << std::endl);
-                maximalDepth = eliminationModelChecker.treatScc(flexibleMatrix, oneStepProbabilities, initialstates, statesToEliminate, forwardTransitions, flexibleBackwardTransitions, false, 0, storm::settings::sparseDtmcEliminationModelCheckerSettings().getMaximalSccSize(), entryStateQueue, missingStateRewards, statePriorities);
-                
-                // If the entry states were to be eliminated last, we need to do so now.
-                STORM_LOG_DEBUG("Eliminating " << entryStateQueue.size() << " entry states as a last step.");
-                if (storm::settings::sparseDtmcEliminationModelCheckerSettings().isEliminateEntryStatesLastSet()) {
-                    for (auto const& state : entryStateQueue) {
-                        eliminationModelChecker.eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
-                    }
-                }
-                subsystem=~statesToEliminate;     
-            }
-            std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << "states." << std::endl;
-            STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " states." << std::endl);
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector<ParametricType>& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix<ParametricType> const& forwardTransitions, boost::optional<std::vector<std::size_t>> const& statePriorities){
-            STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "elimination of states not suported for this type");
-        }
-
-        //OBSOLETE ... 
-        template<>
-        void SparseDtmcRegionModelChecker<storm::RationalFunction, double>::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType>& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities){
-            carl::VariablePool& varPool = carl::VariablePool::getInstance();
-            
-            //first add a state variable for every state in the subsystem, providing that such a variable does not already exist.
-            for (storm::storage::sparse::state_type state : subsystem){
-                if(stateProbVars[state].isZero()){ //variable does not exist yet
-                    storm::Variable stateVar = varPool.getFreshVariable("p_" + std::to_string(state));
-                    std::shared_ptr<carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>> cache(new carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>());
-                    storm::RationalFunction::PolyType stateVarAsPoly(storm::RationalFunction::PolyType::PolyType(stateVar), cache);
-
-                    //each variable is in the interval [0,1]
-                    solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::GEQ, storm::RationalFunction(0));
-                    solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::LEQ, storm::RationalFunction(1));
-                    stateProbVars[state] = stateVarAsPoly;
-                }
-            }
-            
-            //now lets add the actual transitions
-            for (storm::storage::sparse::state_type state : subsystem){
-                storm::RationalFunction reachProbability(oneStepProbabilities[state]);
-                for(auto const& transition : flexibleMatrix.getRow(state)){
-                    reachProbability += transition.getValue() * stateProbVars[transition.getColumn()];
-                }
-                //Todo: depending on the objective (i.e. the formlua) it suffices to use LEQ or GEQ here... maybe this is faster?
-                solver.add(storm::RationalFunction(stateProbVars[state]), storm::CompareRelation::EQ, reachProbability);
-            }
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType>& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities){
-            STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "SMT formulation is not supported for this type");
-        }
-        
-        //DELETEME
-        template<>
-        void SparseDtmcRegionModelChecker<storm::RationalFunction, double>::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType> const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities, ParameterRegion const& region, storm::logic::ComparisonType const& compType){
-            //We are going to build a new (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter 
-                        
-            //todo invent something better to obtain the substitutions.
-            //this only works as long as there is only one parameter per state,
-            // also: check whether the terms are linear/monotone(?)
-            
-            STORM_LOG_WARN("the probability restriction only works on linear terms which is not checked");
-            storm::storage::sparse::state_type const numOfStates=subsystem.getNumberOfSetBits() + 2; //subsystem + target state + sink state
-            storm::models::sparse::StateLabeling stateLabeling(numOfStates);
-            stateLabeling.addLabel("init", storm::storage::BitVector(numOfStates, true));
-            storm::storage::BitVector targetLabel(numOfStates, false);
-            targetLabel.set(numOfStates-2, true);
-            stateLabeling.addLabel("target", std::move(targetLabel));
-            storm::storage::BitVector sinkLabel(numOfStates, false);
-            sinkLabel.set(numOfStates-1, true);
-            stateLabeling.addLabel("sink", std::move(sinkLabel));
-
-            std::pair<storm::storage::SparseMatrix<double>,std::vector<boost::container::flat_set<uint_fast64_t>>> instantiation = instantiateFlexibleMatrix(flexibleMatrix, region.getVerticesOfRegion(region.getVariables()), subsystem, true, oneStepProbabilities, true);
-            boost::optional<std::vector<double>> noStateRewards;
-            boost::optional<storm::storage::SparseMatrix<double>> noTransitionRewards;            
-            storm::models::sparse::Mdp<double> mdp(instantiation.first, std::move(stateLabeling),noStateRewards,noTransitionRewards,instantiation.second);
-                        
-            //we need the correct optimalityType for model checking as well as the correct relation for smt solving
-            storm::logic::OptimalityType opType;
-            storm::CompareRelation boundRelation;
-            switch (compType){
-                case storm::logic::ComparisonType::Greater:
-                    opType=storm::logic::OptimalityType::Minimize;
-                    boundRelation=storm::CompareRelation::GEQ;
-                    break;
-                case storm::logic::ComparisonType::GreaterEqual:
-                    opType=storm::logic::OptimalityType::Minimize;
-                    boundRelation=storm::CompareRelation::GEQ;
-                    break;
-                case storm::logic::ComparisonType::Less:
-                    opType=storm::logic::OptimalityType::Maximize;
-                    boundRelation=storm::CompareRelation::LEQ;
-                    break;
-                case storm::logic::ComparisonType::LessEqual:
-                    opType=storm::logic::OptimalityType::Maximize;
-                    boundRelation=storm::CompareRelation::LEQ;
-                    break;
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-            }
-            
-            //perform model checking on the mdp
-            std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
-            storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr);
-            storm::modelchecker::SparseMdpPrctlModelChecker<double> modelChecker(mdp);
-            std::unique_ptr<CheckResult> resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType);
-            std::vector<double> resultVector = resultPtr->asExplicitQuantitativeCheckResult<double>().getValueVector();            
-            
-            //formulate constraints for the solver
-            uint_fast64_t boundDenominator = 1.0/storm::settings::generalSettings().getPrecision(); //we need to approx. the obtained bounds as rational numbers
-            storm::storage::sparse::state_type subsystemState=0; //the subsystem uses other state indices
-            for(storm::storage::sparse::state_type state : subsystem){
-                uint_fast64_t boundNumerator = resultVector[subsystemState]*boundDenominator;
-                storm::RationalFunction bound(boundNumerator);
-                bound = bound/boundDenominator;
-                //Todo: non-exact values might be problematic here...
-                solver.add(storm::RationalFunction(stateProbVars[state]), boundRelation, bound);
-                ++subsystemState;
-            }
-        }
-        
-        template<typename ParametricType, typename ConstantType>
-        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType> const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities, ParameterRegion const& region, storm::logic::ComparisonType const& compType){
-            STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "restricting Probability Variables is not supported for this type");
-        }
-        
-        //DELETEME
-        template<>
-        bool SparseDtmcRegionModelChecker<storm::RationalFunction, double>::checkRegionOld(storm::logic::Formula const& formula, std::vector<ParameterRegion> parameterRegions){
-            //Note: this is an 'experimental' implementation
-            
-            std::chrono::high_resolution_clock::time_point timeStart = std::chrono::high_resolution_clock::now();
-                
-            //Start with some preprocessing (inspired by computeUntilProbabilities...)
-            //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ])
-            //get the (sub)formulae and the vector of target states
-            STORM_LOG_THROW(formula.isStateFormula(), storm::exceptions::IllegalArgumentException, "expected a stateFormula");
-            STORM_LOG_THROW(formula.asStateFormula().isProbabilityOperatorFormula(), storm::exceptions::IllegalArgumentException, "expected a probabilityOperatorFormula");
-            storm::logic::ProbabilityOperatorFormula const& probOpForm=formula.asStateFormula().asProbabilityOperatorFormula();
-            STORM_LOG_THROW(probOpForm.hasBound(), storm::exceptions::IllegalArgumentException, "The formula has no bound");
-            STORM_LOG_THROW(probOpForm.getSubformula().asPathFormula().isEventuallyFormula(), storm::exceptions::IllegalArgumentException, "expected an eventually subformula");
-            storm::logic::EventuallyFormula const& eventuallyFormula = probOpForm.getSubformula().asPathFormula().asEventuallyFormula();
-            std::unique_ptr<CheckResult> targetStatesResultPtr = eliminationModelChecker.check(eventuallyFormula.getSubformula());
-            storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector();
-            // Do some sanity checks to establish some required properties.
-            STORM_LOG_THROW(model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state.");
-            storm::storage::sparse::state_type initialState = *model.getInitialStates().begin();
-            // Then, compute the subset of states that has a probability of 0 or 1, respectively.
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(model, storm::storage::BitVector(model.getNumberOfStates(),true), targetStates);
-            storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
-            storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
-            storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
-            // If the initial state is known to have either probability 0 or 1, we can directly return the result.
-            if (model.getInitialStates().isDisjointFrom(maybeStates)) {
-                STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step.");
-                double res= statesWithProbability0.get(*model.getInitialStates().begin()) ? 0.0 : 1.0;
-                switch (probOpForm.getComparisonType()){
-                    case storm::logic::ComparisonType::Greater:
-                        return (res > probOpForm.getBound());
-                    case storm::logic::ComparisonType::GreaterEqual:
-                        return (res >= probOpForm.getBound());
-                    case storm::logic::ComparisonType::Less:
-                        return (res < probOpForm.getBound());
-                    case storm::logic::ComparisonType::LessEqual:
-                        return (res <= probOpForm.getBound());
-                    default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-                }
-            }
-            // Determine the set of states that is reachable from the initial state without jumping over a target state.
-            storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, statesWithProbability1);
-            // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
-            maybeStates &= reachableStates;
-            // Create a vector for the probabilities to go to a state with probability 1 in one step.
-            std::vector<storm::RationalFunction> oneStepProbabilities = model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
-            // Determine the set of initial states of the sub-model.
-            storm::storage::BitVector newInitialStates = model.getInitialStates() % maybeStates;
-            // We then build the submatrix that only has the transitions of the maybe states.
-            storm::storage::SparseMatrix<storm::RationalFunction> submatrix = model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
-            storm::storage::SparseMatrix<storm::RationalFunction> submatrixTransposed = submatrix.transpose();
-            
-            std::vector<std::size_t> statePriorities = eliminationModelChecker.getStatePriorities(submatrix, submatrixTransposed, newInitialStates, oneStepProbabilities);
-            
-            // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
-            FlexibleMatrix flexibleMatrix = eliminationModelChecker.getFlexibleSparseMatrix(submatrix);
-            FlexibleMatrix flexibleBackwardTransitions = eliminationModelChecker.getFlexibleSparseMatrix(submatrixTransposed, true);
-            
-            std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now();
-            
-           // Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated.
-            storm::storage::BitVector subsystem = storm::storage::BitVector(submatrix.getRowCount(), true);
-            eliminateStates(subsystem, flexibleMatrix, oneStepProbabilities, flexibleBackwardTransitions, newInitialStates, submatrix, statePriorities);
-            
-            std::chrono::high_resolution_clock::time_point timeStateElemEnd = std::chrono::high_resolution_clock::now();
-            
-            // SMT formulation of resulting pdtmc
-            storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions
-            storm::solver::Smt2SmtSolver solver(manager, true);
-            // we will introduce a variable for every state which encodes the probability to reach a target state from this state.
-            // we will store them as polynomials to easily use operations with rational functions
-            std::vector<storm::RationalFunction::PolyType> stateProbVars(subsystem.size(), storm::RationalFunction::PolyType(0));
-            // todo maybe introduce the parameters already at this point?
-            formulateModelWithSMT(solver, stateProbVars, subsystem, flexibleMatrix, oneStepProbabilities);
-            
-            //the property should be satisfied in the initial state for all parameters.
-            //this is equivalent to:
-            //the negation of the property should not be satisfied for some parameter valuation.
-            //Hence, we flip the comparison relation and later check whether all the constraints are unsat.
-            storm::CompareRelation propertyCompRel;
-            switch (probOpForm.getComparisonType()){
-                case storm::logic::ComparisonType::Greater:
-                    propertyCompRel=storm::CompareRelation::LEQ;
-                    break;
-                case storm::logic::ComparisonType::GreaterEqual:
-                    propertyCompRel=storm::CompareRelation::LT;
-                    break;
-                case storm::logic::ComparisonType::Less:
-                    propertyCompRel=storm::CompareRelation::GEQ;
-                    break;
-                case storm::logic::ComparisonType::LessEqual:
-                    propertyCompRel=storm::CompareRelation::GT;
-                    break;
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
-            }
-            uint_fast64_t thresholdDenominator = 1.0/storm::settings::generalSettings().getPrecision();
-            uint_fast64_t thresholdNumerator = probOpForm.getBound()*thresholdDenominator;
-            storm::RationalFunction threshold(thresholdNumerator);
-            threshold = threshold / thresholdDenominator;
-            solver.add(storm::RationalFunction(stateProbVars[*newInitialStates.begin()]), propertyCompRel, threshold);
-            
-            //the bounds for the parameters
-            solver.push();
-            //STORM_LOG_THROW(parameterRegions.size()==1, storm::exceptions::NotImplementedException, "multiple regions not yet implemented");
-            ParameterRegion region=parameterRegions[0];
-            for(auto variable : region.getVariables()){
-                storm::RawPolynomial lB(variable);
-                lB -= region.getLowerBound(variable);
-                solver.add(carl::Constraint<storm::RawPolynomial>(lB,storm::CompareRelation::GEQ));
-                storm::RawPolynomial uB(variable);
-                uB -= region.getUpperBound(variable);
-                solver.add(carl::Constraint<storm::RawPolynomial>(uB,storm::CompareRelation::LEQ));
-            }
-            
-            std::chrono::high_resolution_clock::time_point timeSmtFormulationEnd = std::chrono::high_resolution_clock::now();
-            
-            // find further restriction on probabilities
-            restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions[0], storm::logic::ComparisonType::Less); //probOpForm.getComparisonType());
-            restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions[0], storm::logic::ComparisonType::Greater);
-            
-            std::chrono::high_resolution_clock::time_point timeRestrictingEnd = std::chrono::high_resolution_clock::now();
-            
-            std::cout << "start solving ..." << std::endl;
-            bool result;
-                switch (solver.check()){
-                case storm::solver::SmtSolver::CheckResult::Sat:
-                    std::cout << "sat!" << std::endl;
-                    result=false;
-                    break;
-                case storm::solver::SmtSolver::CheckResult::Unsat:
-                    std::cout << "unsat!" << std::endl;
-                    result=true;
-                    break;
-                case storm::solver::SmtSolver::CheckResult::Unknown:
-                    std::cout << "unknown!" << std::endl;
-                    STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (Check-result: Unknown)")
-                    result=false;
-                    break;
-                default:
-                    STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (unexpected check-result)")
-                    result=false;
-            }
-            
-            std::chrono::high_resolution_clock::time_point timeSolvingEnd = std::chrono::high_resolution_clock::now();    
-                
-            std::chrono::high_resolution_clock::duration timePreprocessing = timePreprocessingEnd - timeStart;
-            std::chrono::high_resolution_clock::duration timeStateElem = timeStateElemEnd - timePreprocessingEnd;
-            std::chrono::high_resolution_clock::duration timeSmtFormulation = timeSmtFormulationEnd - timeStateElemEnd;
-            std::chrono::high_resolution_clock::duration timeRestricting = timeRestrictingEnd - timeSmtFormulationEnd;
-            std::chrono::high_resolution_clock::duration timeSolving = timeSolvingEnd- timeRestrictingEnd;
-            std::chrono::high_resolution_clock::duration timeOverall = timeSolvingEnd - timeStart;
-            std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timePreprocessing);
-            std::chrono::milliseconds timeStateElemInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeStateElem);
-            std::chrono::milliseconds timeSmtFormulationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeSmtFormulation);
-            std::chrono::milliseconds timeRestrictingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeRestricting);
-            std::chrono::milliseconds timeSolvingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeSolving);
-            std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeOverall);
-            STORM_PRINT_AND_LOG(std::endl << "required time: " << timeOverallInMilliseconds.count() << "ms. Time Breakdown:" << std::endl);
-            STORM_PRINT_AND_LOG("    * " << timePreprocessingInMilliseconds.count() << "ms for Preprocessing" << std::endl);
-            STORM_PRINT_AND_LOG("    * " << timeStateElemInMilliseconds.count() << "ms for StateElemination" << std::endl);
-            STORM_PRINT_AND_LOG("    * " << timeSmtFormulationInMilliseconds.count() << "ms for SmtFormulation" << std::endl);
-            STORM_PRINT_AND_LOG("    * " << timeRestrictingInMilliseconds.count() << "ms for Restricting" << std::endl);
-            STORM_PRINT_AND_LOG("    * " << timeSolvingInMilliseconds.count() << "ms for Solving" << std::endl);
-
-            return result;
-        }
-
-        template<typename ParametricType, typename ConstantType>
-        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkRegionOld(storm::logic::Formula const& formula, std::vector<ParameterRegion> parameterRegions){
-            STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Region check is not supported for this type");
-        }
-        
-        #endif
-
-
-        
-        
-#ifdef STORM_HAVE_CARL
-        template class SparseDtmcRegionModelChecker<storm::RationalFunction, double>;
-#endif
-        //note: for other template instantiations, add a rule for the typedefs of VariableType and CoefficientType
-        
-    } // namespace modelchecker
-} // namespace storm
diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h
deleted file mode 100644
index 7f58754f3..000000000
--- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h
+++ /dev/null
@@ -1,386 +0,0 @@
-#ifndef STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_
-#define STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_
-
-#include <type_traits>
-
-#include "src/storage/sparse/StateType.h"
-#include "src/models/sparse/Dtmc.h"
-#include "src/models/sparse/Mdp.h"
-#include "src/utility/constants.h"
-#include "utility/regions.h"
-#include "src/solver/Smt2SmtSolver.h"
-#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
-
-
-
-namespace storm {
-    namespace modelchecker {
-        
-        template<typename ParametricType, typename ConstantType>
-        class SparseDtmcRegionModelChecker {
-        public:
-            
-           
-            
-            //The type of variables and bounds depends on the template arguments
-            typedef typename std::conditional<(std::is_same<ParametricType,storm::RationalFunction>::value), storm::Variable,std::nullptr_t>::type VariableType;
-            typedef typename std::conditional<(std::is_same<ParametricType,storm::RationalFunction>::value), storm::Coefficient,std::nullptr_t>::type CoefficientType;
-            
-            /*!
-             * The possible results for a single region
-             */
-            enum class RegionCheckResult { 
-                UNKNOWN, /*!< the result is unknown */
-                EXISTSSAT, /*!< the formula is satisfied for at least one parameter evaluation that lies in the given region */
-                EXISTSVIOLATED, /*!< the formula is violated for at least one parameter evaluation that lies in the given region */
-                EXISTSBOTH, /*!< the formula is satisfied for some parameters but also violated for others */
-                ALLSAT, /*!< the formula is satisfied for all parameters in the given region */
-                ALLVIOLATED /*!< the formula is violated for all parameters in the given region */
-            };
-            
-            class ApproximationModel;
-            class SamplingModel;
-            
-            class ParameterRegion{
-            public:
-
-                ParameterRegion(std::map<VariableType, CoefficientType> lowerBounds, std::map<VariableType, CoefficientType> upperBounds);
-                
-                
-                
-                std::set<VariableType> getVariables() const;
-                CoefficientType const& getLowerBound(VariableType const& variable) const;
-                CoefficientType const& getUpperBound(VariableType const& variable) const;
-                const std::map<VariableType, CoefficientType> getUpperBounds() const;
-                const std::map<VariableType, CoefficientType> getLowerBounds() const;
-                
-                /*
-                 * Returns a vector of all possible combinations of lower and upper bounds of the given variables.
-                 * The first entry of the returned vector will map every variable to its lower bound
-                 * The second entry will map every variable to its lower bound, except the first one (i.e. *consVariables.begin())
-                 * ...
-                 * The last entry will map every variable to its upper bound
-                 * 
-                 * If the given set of variables is empty, the returned vector will contain an empty map
-                 */
-                std::vector<std::map<VariableType, CoefficientType>> getVerticesOfRegion(std::set<VariableType> const& consideredVariables) const;
-                
-                //returns the currently set check result as a string
-                std::string checkResultToString() const;
-                
-                //returns the region as string in the format 0.3<=p<=0.4,0.2<=q<=0.5;
-                std::string toString() const;
-                
-                void setCheckResult(RegionCheckResult checkResult);
-                RegionCheckResult getCheckResult() const;
-                
-                /*!
-                 * Sets a point in the region for which the considered property is not satisfied. 
-                 */
-                void setViolatedPoint(std::map<VariableType, CoefficientType> const& violatedPoint);
-                
-                /*!
-                 * Retrieves a point in the region for which is considered property is not satisfied.
-                 * If such a point is not known, the returned map is empty.
-                 */
-                std::map<VariableType, CoefficientType> getViolatedPoint() const;
-                
-                
-                /*!
-                 * Sets a point in the region for which the considered property is satisfied. 
-                 */
-                void setSatPoint(std::map<VariableType, CoefficientType> const& satPoint);
-                
-                /*!
-                 * Retrieves a point in the region for which is considered property is satisfied.
-                 * If such a point is not known, the returned map is empty.
-                 */
-                std::map<VariableType, CoefficientType> getSatPoint() const;
-                
-            private:
-                
-                std::map<VariableType, CoefficientType> const lowerBounds;
-                std::map<VariableType, CoefficientType> const upperBounds;
-                std::set<VariableType> variables;
-                RegionCheckResult checkResult;
-                std::map<VariableType, CoefficientType> satPoint;
-                std::map<VariableType, CoefficientType> violatedPoint;
-                
-                
-            };
-            
-            explicit SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc<ParametricType> const& model);
-
-            /*!
-             * Checks if the given formula can be handled by this
-             * @param formula the formula to be checked
-             */
-            bool canHandle(storm::logic::Formula const& formula) const;
-            
-            /*!
-             * Specifies the considered formula.
-             * A few preprocessing steps are performed.
-             * If another formula has been specified before, all 'context' regarding the old formula is lost.
-             * 
-             * @param formula the formula to be considered.
-             */
-            void specifyFormula(storm::logic::Formula const& formula);
-
-            /*!
-             * Checks whether the given formula holds for all parameters that lie in the given region.
-             * Sets the region checkresult accordingly. Moreover, region.satPoint and/or an region.violatedPoint will be set.
-             * 
-             * @note A formula has to be specified first.
-             * 
-             * @param region The considered region
-             * 
-             */
-            void checkRegion(ParameterRegion& region);
-            
-            /*!
-             * Checks for every given region whether the specified formula holds for all parameters that lie in that region.
-             * Sets the region checkresult accordingly. Moreover, region.satPoint and/or an region.violatedPoint will be set.
-             * 
-             * @note A formula has to be specified first.
-             * 
-             * @param region The considered region
-             */
-            void checkRegions(std::vector<ParameterRegion>& regions);
-            
-            /*!
-             * Checks whether the given formula holds for all possible parameters that satisfy the given parameter regions
-             * ParameterRegions should contain all parameters.
-             */
-            bool checkRegionOld(storm::logic::Formula const& formula, std::vector<ParameterRegion> parameterRegions);
-            
-            /*!
-             * Prints statistical information (mostly running times) to the given stream.
-             */
-            void printStatisticsToStream(std::ostream& outstream);
-            
-            
-            
-        private:
-            
-            typedef typename storm::modelchecker::SparseDtmcEliminationModelChecker<ParametricType>::FlexibleSparseMatrix FlexibleMatrix;
-            
-            /*!
-             * Represents the current state of this
-             */
-           // enum class RegionCheckerState{
-         //       NOFORMULA, /*!< there is no formula */
-         //       INITIALIZED, /*!< a formula has been specified. Ready to get regions*/
-        //    };
-            
-   #ifdef STORM_HAVE_CARL
-                /*!
-                 * Instantiates the matrix, i.e., evaluate the occurring functions according to the given substitutions of the variables.
-                 * One row of the given flexible matrix will be one rowgroup in the returned matrix, consisting of one row for every substitution.
-                 * The returned matrix can be seen as the transition matrix of an MDP with the action labeling given by the returned vector of sets.
-                 * Only the rows selected by the given filter are considered. (filter should have size==matrix.getNumberOfRows())
-                 * An exception is thrown if there is a transition from a selected state to an unselected state
-                 * If one step probabilities are given, a new state is added which can be considered as target state.
-                 * The "missing" probability can be redirected to a sink state
-                 * By convention, the target state will have index filter.getNumberOfSetBits() and the sink state will be the state with the highest index (so right after the target state)
-                 * 
-                 * @param matrix the considered flexible matrix
-                 * @param substitutions A list of mappings, each assigning a constant value to every variable
-                 * @param filter selects the rows of this flexibleMatrix, that will be considered
-                 * @param addSinkState adds a state with a self loop to which the "missing" probability will lead
-                 * @param oneStepProbabilities if given, a new state is added to which there are transitions for all non-zero entries in this vector
-                 * @param addSelfLoops if set, zero valued selfloops will be added in every row
-                 * 
-                 * @return A matrix with constant (double) entries and a choice labeling
-                 */
-                std::pair<storm::storage::SparseMatrix<double>,std::vector<boost::container::flat_set<uint_fast64_t>>> instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector<std::map<storm::Variable, storm::RationalFunction::CoeffType>> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState=true, std::vector<ParametricType> const& oneStepProbabilities=std::vector<ParametricType>(), bool addSelfLoops=true) const;
-                //todo add const keyword
-            
-            //eliminates some of the states according to different strategies.
-            void eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector<ParametricType>& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::SparseMatrix<ParametricType> const& forwardTransitions, boost::optional<std::vector<std::size_t>> const& statePriorities = {});
-            
-            void formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType>& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities);
-            
-            void restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector<storm::RationalFunction::PolyType> const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector<storm::RationalFunction> const& oneStepProbabilities, ParameterRegion const& region, storm::logic::ComparisonType const& compTypeOfProperty);
-#endif         
-            template <typename ValueType>
-            bool valueIsInBoundOfFormula(ValueType value);
-            
-            /*!
-             * eliminates all states for which the outgoing transitions are constant.
-             * Also checks whether the non constant functions are linear
-             */
-            void computeSimplifiedModel(storm::storage::BitVector const& targetStates);
-            
-            /*!
-             * Initializes a Sample-Model that can be used to get the probability result for a certain parameter evaluation.
-             * Initializes an Approximation-Model that can be used to approximate the reachability probabilities.
-             */
-            void initializeSampleAndApproxModel();
-            
-            /*!
-             * Initializes a DTMC that can be used to get the probability result for a certain parameter evaluation.
-             * To quickly insert different evaluations, we provide a dtmc with some dummy values in the transition entries.
-             * Furthermore, another matrix is given that is used to connect matrix entries and transition functions.
-             * In a similar way an MDP is initialized. The Mdp can be used to approximate the probabilities.
-             * In addition to matrix entries and transition functions, we specify which variables and which bounds needs to be substituted
-             * 
-             * 
-             * @param sampleDtmc This dtmc can represent different instantiations of the considered pDtmc
-             * @param sampleDtmcMapping Connects the entries of the sampleDtmc Matrix with the corresponding transitionFunctions
-             * @param approxMdp This mdp will later be used to approximate the reachability probabilities
-             * @param approxMdpMapping Connects the entries of the approxMdp Matrix with the corresponding transitionFunctions and a substitution (given as an index in the approxMdpSubstitutionsVector)
-             * @param approxMdpSubstitutions contains the substitutions of the parameters
-             * 
-             * @param subsys the states of the flexTransitions that are still part of the pDTMC (i.e. that have not been eliminated)
-             * @param flexTransitions the transitions of the pDTMC 
-             * @param oneStepProbs the probabilities to move to a target state
-             * @param initState the initial state of the pDtmc
-            void initializeSampleDtmcAndApproxMdp2(
-                    std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> & sampleDtmc,
-                    std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>> & sampleDtmcMapping,
-                    std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> & approxMdp,
-                    std::vector<std::tuple<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&, std::size_t>> & approxMdpMapping,
-                    std::vector<std::map<VariableType, TypeOfBound>> & approxMdpSubstitutions,
-                    storm::storage::BitVector const& subsys,
-                    storm::storage::SparseMatrix<ParametricType> const& transitions,
-                    std::vector<ParametricType> const& oneStepProbs,
-                    storm::storage::sparse::state_type const& initState
-            );
-             */
-            
-            ParametricType getReachProbFunction();
-            
-            
-            //initializes the given solver which can later be used to give an exact result regarding the whole model.
-            void initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, ParametricType const& reachProbFunction, storm::logic::ProbabilityOperatorFormula const& formula);
-            
-            /*!
-             * Checks the value of the function at some sampling points within the given region.
-             * May set the satPoint and violatedPoint of the regions if they are not yet specified and such points are found
-             * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
-             * 
-             * @return true if an violated point as well as a sat point has been found during the process
-             */
-            bool checkSamplePoints(ParameterRegion& region);
-            
-            /*!
-             * Checks the value of the function at the given sampling point.
-             * May set the satPoint and violatedPoint of the regions if thy are not yet specified and such point is given.
-             * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
-             * 
-             * @param viaReachProbFunction if set, the sampling will be done via the reachProbFunction.
-             *                             Otherwise, the model will be instantiated and checked
-             * 
-             * @return true if an violated point as well as a sat point has been found, i.e., the check result is changed to EXISTSOTH
-             */
-            bool checkPoint(ParameterRegion& region, std::map<VariableType, CoefficientType>const& point, bool viaReachProbFunction=false);
-            
-            /*!
-             * Builds an MDP that is used to compute bounds on the maximal/minimal reachability probability,
-             * If this approximation already yields that the property is satisfied/violated in the whole region,
-             * true is returned and the regioncheckresult of the region is changed accordingly.
-             * The computed bounds are stored in the given vectors.
-             * However, if only the lowerBounds (or upperBounds) need to be computed in order to prove ALLSAT or ALLVIOLATED, 
-             * the upperBounds (or lowerBounds) vector remains untouched.
-             */
-            bool checkApproximativeProbabilities(ParameterRegion& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds); 
-            
-            
-            /*!
-             * Builds the mdp that is used to obtain bounds on the maximal/minimal reachability probability
-             * The result is stored in this->approxMdp
-             */
-            void buildMdpForApproximation(ParameterRegion const& region);
-            
-            /*!
-             * Builds the mdp that is used to obtain bounds on the maximal/minimal reachability probability
-             
-             */
-            storm::models::sparse::Mdp<ConstantType> buildMdpForApproximation2(ParameterRegion const& region);
-            
-            /*!
-             * Starts the SMTSolver to get the result.
-             * The current regioncheckresult of the region should be EXISTSSAT or EXISTVIOLATED.
-             * Otherwise, a sampingPoint will be computed.
-             * True is returned iff the solver was successful (i.e., it returned sat or unsat)
-             * A Sat- or Violated point is set, if the solver has found one.
-             * The region checkResult of the given region is changed accordingly.
-             */
-            bool checkFullSmt(ParameterRegion& region); 
-            
-            
-            // The model this model checker is supposed to analyze.
-            storm::models::sparse::Dtmc<ParametricType> const& model;
-
-            //classes that provide auxilliary functions
-            // Instance of an elimination model checker to access its functions
-            storm::modelchecker::SparseDtmcEliminationModelChecker<ParametricType> eliminationModelChecker;
-            // comparators that can be used to compare constants.
-            storm::utility::ConstantsComparator<ParametricType> parametricTypeComparator;
-            storm::utility::ConstantsComparator<ConstantType> constantTypeComparator;
-            
-            
-            std::shared_ptr<storm::solver::Smt2SmtSolver> smtSolver;
-            
-
-            //the following members depend on the currently specified formula:
-            
-            //the currently specified formula
-            std::unique_ptr<storm::logic::ProbabilityOperatorFormula> probabilityOperatorFormula;
-            
-            // the original model after states with constant transitions have been eliminated
-            std::shared_ptr<storm::models::sparse::Dtmc<ParametricType>> simplifiedModel;
-            
-            // a flag that is true if there are only linear functions at transitions of the model
-            bool hasOnlyLinearFunctions;
-            
-            // the model that can be instantiated to check the value at a certain point
-            std::shared_ptr<SamplingModel> samplingModel;
-            // the model that  is used to approximate the probability values
-            std::shared_ptr<ApproximationModel> approximationModel;
-            
-            /*
-            // the dtmc that can be instantiated to check the value at a certain point
-            std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> sampleDtmc;
-            // a vector that links entries of the dtmc matrix with the corresponding transition functions (for fast instantiation)
-            std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>> sampleDtmcMapping;
-            // the Mdp that is used to approximate the probability values
-            std::shared_ptr<storm::models::sparse::Mdp<ConstantType>>  approxMdp;
-            // a vector that links entries of the mdp matrix with the corresponding transition functions and substitutions (given as an index of this->approxMdpSubstitutions vector)
-            std::vector<std::tuple<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&, std::size_t>> approxMdpMapping;
-            // stores the different substitutions of the variables
-            std::vector<std::map<VariableType, TypeOfBound>> approxMdpSubstitutions;
-              */      
-            // The  function for the reachability probability in the initial state 
-            ParametricType reachProbFunction;
-            bool isReachProbFunctionComputed;
-            bool isResultConstant;
-            
-            
-            // runtimes and other information for statistics. 
-            uint_fast64_t numOfCheckedRegions;
-            
-            uint_fast64_t numOfRegionsSolvedThroughSampling;
-            uint_fast64_t numOfRegionsSolvedThroughApproximation;
-            uint_fast64_t numOfRegionsSolvedThroughSubsystemSmt;
-            uint_fast64_t numOfRegionsSolvedThroughFullSmt;
-            
-            uint_fast64_t numOfRegionsExistsBoth;
-            uint_fast64_t numOfRegionsAllSat;
-            uint_fast64_t numOfRegionsAllViolated;
-            
-            std::chrono::high_resolution_clock::duration timePreprocessing;
-            std::chrono::high_resolution_clock::duration timeInitialStateElimination;
-            std::chrono::high_resolution_clock::duration timeComputeReachProbFunction;
-            std::chrono::high_resolution_clock::duration timeCheckRegion;
-            std::chrono::high_resolution_clock::duration timeSampling;
-            std::chrono::high_resolution_clock::duration timeApproximation;
-            std::chrono::high_resolution_clock::duration timeMDPBuild;
-            std::chrono::high_resolution_clock::duration timeSubsystemSmt;
-            std::chrono::high_resolution_clock::duration timeFullSmt;
-        };
-        
-    } // namespace modelchecker
-} // namespace storm
-
-#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_ */
diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp
index 3e157e62b..ae155fb53 100644
--- a/src/modelchecker/region/ApproximationModel.cpp
+++ b/src/modelchecker/region/ApproximationModel.cpp
@@ -6,6 +6,7 @@
  */
 
 #include "src/modelchecker/region/ApproximationModel.h"
+#include "src/modelchecker/region/ParameterRegion.h"
 #include "modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "exceptions/UnexpectedException.h"
@@ -61,7 +62,7 @@ namespace storm {
                     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));
+                            auto pair=distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex));
                             ++numOfNonConstEntries;
                         }
                         matrixBuilder.addNextValue(approxModelRow, entry.getColumn(), dummyEntry);
@@ -92,11 +93,16 @@ namespace storm {
                     auto appEntry = this->model->getTransitionMatrix().getRow(approxModelRow).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<ParametricType, ConstantType>(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[approxModelRow], storm::utility::zero<ConstantType>());
-                            auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple);
+                            auto const tableIt= std::find(evaluationTable.begin(), evaluationTable.end(), searchedTuple);
+                            //auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple);
+                            //auto const& tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple);
+                            if(tableIt==evaluationTable.end()){
+                                std::cout << "did not found tuple in the table: " << parEntry.getValue() << " substitution " << rowSubstitutions[approxModelRow] << " in parametric model at row " << row << " column " << parEntry.getColumn() << " approximation Row " << approxModelRow << std::endl;
+                            }
                             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");
                             mapping.emplace_back(std::make_pair(&(std::get<2>(*tableIt)), appEntry));
                         }
@@ -138,7 +144,7 @@ namespace storm {
             //write entries into evaluation table
             for(auto& tableEntry : this->evaluationTable){
                 std::get<2>(tableEntry)=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
-                        storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(
+                        storm::utility::regions::evaluateFunction(
                                 std::get<0>(tableEntry),
                                 instantiatedSubs[std::get<1>(tableEntry)]
                             )
diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h
index b2da957b9..06bc06155 100644
--- a/src/modelchecker/region/ApproximationModel.h
+++ b/src/modelchecker/region/ApproximationModel.h
@@ -8,7 +8,8 @@
 #ifndef STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H
 #define	STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H
 
-#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
+#include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
+#include "src/models/sparse/Mdp.h"
 
 namespace storm {
     namespace modelchecker {
diff --git a/src/modelchecker/region/ParameterRegion.cpp b/src/modelchecker/region/ParameterRegion.cpp
new file mode 100644
index 000000000..6bc33ccc5
--- /dev/null
+++ b/src/modelchecker/region/ParameterRegion.cpp
@@ -0,0 +1,258 @@
+/* 
+ * File:   ParameterRegion.cpp
+ * Author: tim
+ * 
+ * Created on August 10, 2015, 1:51 PM
+ */
+
+#include "src/modelchecker/region/ParameterRegion.h"
+
+#include "src/utility/regions.h"
+
+#include "src/exceptions/UnexpectedException.h"
+#include "exceptions/InvalidSettingsException.h"
+#include "parser/MappedFile.h"
+
+namespace storm {
+    namespace modelchecker {
+
+        template<typename ParametricType, typename ConstantType>
+        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::ParameterRegion(std::map<VariableType, CoefficientType> lowerBounds, std::map<VariableType, CoefficientType> upperBounds) : lowerBounds(lowerBounds), upperBounds(upperBounds), checkResult(RegionCheckResult::UNKNOWN) {
+            //check whether both mappings map the same variables and precompute the set of variables
+            for (auto const& variableWithBound : lowerBounds) {
+                STORM_LOG_THROW((upperBounds.find(variableWithBound.first) != upperBounds.end()), storm::exceptions::InvalidArgumentException, "Couldn't create region. No upper bound specified for Variable " << variableWithBound.first);
+                this->variables.insert(variableWithBound.first);
+            }
+            for (auto const& variableWithBound : upperBounds) {
+                STORM_LOG_THROW((this->variables.find(variableWithBound.first) != this->variables.end()), storm::exceptions::InvalidArgumentException, "Couldn't create region. No lower bound specified for Variable " << variableWithBound.first);
+            }
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::~ParameterRegion() {
+            //Intentionally left empty
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        std::set<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getVariables() const {
+            return this->variables;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getLowerBound(VariableType const& variable) const {
+            auto const& result = lowerBounds.find(variable);
+            STORM_LOG_THROW(result != lowerBounds.end(), storm::exceptions::IllegalArgumentException, "tried to find a lower bound of a variable that is not specified by this region");
+            return (*result).second;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getUpperBound(VariableType const& variable) const {
+            auto const& result = upperBounds.find(variable);
+            STORM_LOG_THROW(result != upperBounds.end(), storm::exceptions::IllegalArgumentException, "tried to find an upper bound of a variable that is not specified by this region");
+            return (*result).second;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        const std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getUpperBounds() const {
+            return upperBounds;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        const std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getLowerBounds() const {
+            return lowerBounds;
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        std::vector<std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType>> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getVerticesOfRegion(std::set<VariableType> const& consideredVariables) const {
+            std::size_t const numOfVariables = consideredVariables.size();
+            std::size_t const numOfVertices = std::pow(2, numOfVariables);
+            std::vector<std::map<VariableType, CoefficientType >> resultingVector(numOfVertices, std::map<VariableType, CoefficientType>());
+            if (numOfVertices == 1) {
+                //no variables are given, the returned vector should still contain an empty map
+                return resultingVector;
+            }
+
+            for (uint_fast64_t vertexId = 0; vertexId < numOfVertices; ++vertexId) {
+                //interprete vertexId as a bit sequence
+                //the consideredVariables.size() least significant bits of vertex will always represent the next vertex
+                //(00...0 = lower bounds for all variables, 11...1 = upper bounds for all variables)
+                std::size_t variableIndex = 0;
+                for (auto const& variable : consideredVariables) {
+                    if ((vertexId >> variableIndex) % 2 == 0) {
+                        resultingVector[vertexId].insert(std::pair<VariableType, CoefficientType>(variable, getLowerBound(variable)));
+                    } else {
+                        resultingVector[vertexId].insert(std::pair<VariableType, CoefficientType>(variable, getUpperBound(variable)));
+                    }
+                    ++variableIndex;
+                }
+            }
+            return resultingVector;
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::RegionCheckResult SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getCheckResult() const {
+            return checkResult;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setCheckResult(RegionCheckResult checkResult) {
+            //a few sanity checks
+            STORM_LOG_THROW((this->checkResult == RegionCheckResult::UNKNOWN || checkResult != RegionCheckResult::UNKNOWN), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from something known to UNKNOWN ");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSSAT || checkResult != RegionCheckResult::EXISTSVIOLATED), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSSAT to EXISTSVIOLATED");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSSAT || checkResult != RegionCheckResult::ALLVIOLATED), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSSAT to ALLVIOLATED");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSVIOLATED || checkResult != RegionCheckResult::EXISTSSAT), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSVIOLATED to EXISTSSAT");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSVIOLATED || checkResult != RegionCheckResult::ALLSAT), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSVIOLATED to ALLSAT");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSBOTH || checkResult != RegionCheckResult::ALLSAT), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSBOTH to ALLSAT");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::EXISTSBOTH || checkResult != RegionCheckResult::ALLVIOLATED), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from EXISTSBOTH to ALLVIOLATED");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::ALLSAT || checkResult == RegionCheckResult::ALLSAT), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from ALLSAT to something else");
+            STORM_LOG_THROW((this->checkResult != RegionCheckResult::ALLVIOLATED || checkResult == RegionCheckResult::ALLVIOLATED), storm::exceptions::InvalidArgumentException, "Tried to change the check result of a region from ALLVIOLATED to something else");
+            this->checkResult = checkResult;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getViolatedPoint() const {
+            return violatedPoint;
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setViolatedPoint(std::map<VariableType, CoefficientType> const& violatedPoint) {
+            this->violatedPoint = violatedPoint;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        std::map<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType, typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getSatPoint() const {
+            return satPoint;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::setSatPoint(std::map<VariableType, CoefficientType> const& satPoint) {
+            this->satPoint = satPoint;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        std::string SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::checkResultToString() const {
+            switch (this->checkResult) {
+                case RegionCheckResult::UNKNOWN:
+                    return "unknown";
+                case RegionCheckResult::EXISTSSAT:
+                    return "ExistsSat";
+                case RegionCheckResult::EXISTSVIOLATED:
+                    return "ExistsViolated";
+                case RegionCheckResult::EXISTSBOTH:
+                    return "ExistsBoth";
+                case RegionCheckResult::ALLSAT:
+                    return "allSat";
+                case RegionCheckResult::ALLVIOLATED:
+                    return "allViolated";
+            }
+            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not identify check result")
+            return "ERROR";
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        std::string SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::toString() const {
+            std::stringstream regionstringstream;
+            for (auto var : this->getVariables()) {
+                regionstringstream << storm::utility::regions::convertNumber<SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType, double>(this->getLowerBound(var));
+                regionstringstream << "<=";
+                regionstringstream << storm::utility::regions::getVariableName(var);
+                regionstringstream << "<=";
+                regionstringstream << storm::utility::regions::convertNumber<SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType, double>(this->getUpperBound(var));
+                regionstringstream << ",";
+            }
+            std::string regionstring = regionstringstream.str();
+            //the last comma should actually be a semicolon
+            regionstring = regionstring.substr(0, regionstring.length() - 1) + ";";
+            return regionstring;
+        }
+
+        
+        
+              template<typename ParametricType, typename ConstantType>
+            void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::parseParameterBounds(
+                    std::map<VariableType, CoefficientType>& lowerBounds,
+                    std::map<VariableType, CoefficientType>& upperBounds,
+                    std::string const& parameterBoundsString,
+                    double const precision){
+                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
+                
+                std::string::size_type positionOfFirstRelation = parameterBoundsString.find("<=");
+                STORM_LOG_THROW(positionOfFirstRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find  a '<=' after the first number");
+                std::string::size_type positionOfSecondRelation = parameterBoundsString.find("<=", positionOfFirstRelation+2);
+                STORM_LOG_THROW(positionOfSecondRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find  a '<=' after the parameter");
+                
+                std::string parameter=parameterBoundsString.substr(positionOfFirstRelation+2,positionOfSecondRelation-(positionOfFirstRelation+2));
+                //removes all whitespaces from the parameter string:
+                parameter.erase(std::remove_if(parameter.begin(), parameter.end(), ::isspace), parameter.end());
+                STORM_LOG_THROW(parameter.length()>0, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find a parameter");
+                double lowerBound, upperBound;
+                try{
+                    lowerBound=std::stod(parameterBoundsString.substr(0,positionOfFirstRelation));
+                    upperBound=std::stod(parameterBoundsString.substr(positionOfSecondRelation+2));
+                }
+                catch (std::exception const& exception) {
+                    STORM_LOG_ERROR("Failed to parse the region: " << parameterBoundsString << ". The correct format for regions is lowerBound<=parameter<=upperbound");
+                    throw exception;
+                }
+                
+                VariableType var = storm::utility::regions::getVariableFromString<VariableType>(parameter);
+                CoefficientType lb = storm::utility::regions::convertNumber<double, CoefficientType>(lowerBound, true, actualPrecision);
+                STORM_LOG_WARN_COND((lowerBound==storm::utility::regions::convertNumber<CoefficientType, double>(lb, true, actualPrecision)), "The lower bound of '"<< parameterBoundsString << "' could not be parsed accurately. Increase precision?");
+                CoefficientType ub = storm::utility::regions::convertNumber<double, CoefficientType>(upperBound, false, actualPrecision);
+                STORM_LOG_WARN_COND((upperBound==storm::utility::regions::convertNumber<CoefficientType, double>(ub, true, actualPrecision)), "The upper bound of '"<< parameterBoundsString << "' could not be parsed accurately. Increase precision?");
+                lowerBounds.emplace(std::make_pair(var, lb));  
+                upperBounds.emplace(std::make_pair(var, ub));
+               // std::cout << "parsed bounds " << parameterBoundsString << ": lb=" << lowerBound << " ub=" << upperBound << " param='" << parameter << "' precision=" << actualPrecision << std::endl;
+            }
+            
+            template<typename ParametricType, typename ConstantType>
+            typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::parseRegion(std::string const& regionString, double precision){
+                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
+                std::map<VariableType, CoefficientType> lowerBounds;
+                std::map<VariableType, CoefficientType> upperBounds;
+                std::vector<std::string> parameterBounds;
+                boost::split(parameterBounds, regionString, boost::is_any_of(","));
+                for(auto const& parameterBound : parameterBounds){
+                    parseParameterBounds(lowerBounds, upperBounds, parameterBound, actualPrecision);
+                }
+                return ParameterRegion(lowerBounds, upperBounds);
+            }
+            
+            template<typename ParametricType, typename ConstantType>
+            std::vector<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::parseMultipleRegions(std::string const& regionsString, double precision){
+                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
+                std::vector<ParameterRegion> result;
+                std::vector<std::string> regionsStrVec;
+                boost::split(regionsStrVec, regionsString, boost::is_any_of(";"));
+                for(auto const& regionStr : regionsStrVec){
+                    if(!std::all_of(regionStr.begin(),regionStr.end(), ::isspace)){ //skip this string if it only consists of space
+                    result.emplace_back(parseRegion(regionStr, actualPrecision));
+                    }
+                }
+                return result;
+            }
+            
+            template<typename ParametricType, typename ConstantType>
+            std::vector<typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion::getRegionsFromSettings(double precision){
+                STORM_LOG_THROW(storm::settings::regionSettings().isRegionsSet() || storm::settings::regionSettings().isRegionFileSet(), storm::exceptions::InvalidSettingsException, "Tried to obtain regions from the settings but no regions are specified.");
+                STORM_LOG_THROW(!(storm::settings::regionSettings().isRegionsSet() && storm::settings::regionSettings().isRegionFileSet()), storm::exceptions::InvalidSettingsException, "Regions are specified via file AND cmd line. Only one option is allowed.");
+                
+                std::string regionsString;
+                if(storm::settings::regionSettings().isRegionsSet()){
+                    regionsString = storm::settings::regionSettings().getRegionsFromCmdLine();
+                }
+                else{
+                    //if we reach this point we can assume that the region is given as a file.
+                    STORM_LOG_THROW(storm::parser::MappedFile::fileExistsAndIsReadable(storm::settings::regionSettings().getRegionFilePath().c_str()), storm::exceptions::InvalidSettingsException, "The path to the file in which the regions are specified is not valid.");
+                    storm::parser::MappedFile mf(storm::settings::regionSettings().getRegionFilePath().c_str());
+                    regionsString = std::string(mf.getData(),mf.getDataSize());
+                }
+                return parseMultipleRegions(regionsString,precision);
+            }
+#ifdef STORM_HAVE_CARL
+        template class SparseDtmcRegionModelChecker<storm::RationalFunction, double>::ParameterRegion;
+#endif
+
+    }
+}
+
diff --git a/src/modelchecker/region/ParameterRegion.h b/src/modelchecker/region/ParameterRegion.h
new file mode 100644
index 000000000..982b67c55
--- /dev/null
+++ b/src/modelchecker/region/ParameterRegion.h
@@ -0,0 +1,136 @@
+/* 
+ * File:   ParameterRegion.h
+ * Author: tim
+ *
+ * Created on August 10, 2015, 1:51 PM
+ */
+
+#ifndef STORM_MODELCHECKER_REGION_PARAMETERREGION_H
+#define	STORM_MODELCHECKER_REGION_PARAMETERREGION_H
+
+#include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
+
+namespace storm {
+    namespace modelchecker{
+
+        template<typename ParametricType, typename ConstantType>
+        class SparseDtmcRegionModelChecker;
+        
+        template<typename ParametricType, typename ConstantType>
+        class SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ParameterRegion{
+        public:
+            typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType VariableType;
+            typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType CoefficientType;
+            
+            ParameterRegion(std::map<VariableType, CoefficientType> lowerBounds, std::map<VariableType, CoefficientType> upperBounds);
+            virtual ~ParameterRegion();
+                
+            std::set<VariableType> getVariables() const;
+            CoefficientType const& getLowerBound(VariableType const& variable) const;
+            CoefficientType const& getUpperBound(VariableType const& variable) const;
+            const std::map<VariableType, CoefficientType> getUpperBounds() const;
+            const std::map<VariableType, CoefficientType> getLowerBounds() const;
+                
+            /*
+             * Returns a vector of all possible combinations of lower and upper bounds of the given variables.
+             * The first entry of the returned vector will map every variable to its lower bound
+             * The second entry will map every variable to its lower bound, except the first one (i.e. *getVariables.begin())
+             * ...
+             * The last entry will map every variable to its upper bound
+             * 
+             * If the given set of variables is empty, the returned vector will contain an empty map
+             */
+            std::vector<std::map<VariableType, CoefficientType>> getVerticesOfRegion(std::set<VariableType> const& consideredVariables) const;
+         
+            RegionCheckResult getCheckResult() const;
+            void setCheckResult(RegionCheckResult checkResult);
+ 
+            /*!
+             * Retrieves a point in the region for which is considered property is not satisfied.
+             * If such a point is not known, the returned map is empty.
+             */
+            std::map<VariableType, CoefficientType> getViolatedPoint() const;
+
+            /*!
+             * Sets a point in the region for which the considered property is not satisfied. 
+             */
+            void setViolatedPoint(std::map<VariableType, CoefficientType> const& violatedPoint);
+                
+            /*!
+             * Retrieves a point in the region for which is considered property is satisfied.
+             * If such a point is not known, the returned map is empty.
+             */
+            std::map<VariableType, CoefficientType> getSatPoint() const;
+                
+            /*!
+             * Sets a point in the region for which the considered property is satisfied. 
+             */
+            void setSatPoint(std::map<VariableType, CoefficientType> const& satPoint);
+            
+            //returns the currently set check result as a string
+            std::string checkResultToString() const;
+     
+            //returns the region as string in the format 0.3<=p<=0.4,0.2<=q<=0.5;
+            std::string toString() const;
+
+            /*
+             * Can be used to parse a single parameter with its bounds from a string of the form "0.3<=p<=0.5".
+             * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
+             * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
+             * If no precision is given, the one from the settings is used.
+             * The results will be inserted in the given maps
+             * 
+             */
+            static void parseParameterBounds( 
+                    std::map<VariableType, CoefficientType>& lowerBounds,
+                    std::map<VariableType, CoefficientType>& upperBounds,
+                    std::string const& parameterBoundsString,
+                    double const precision=0.0
+            );
+
+            /*
+             * Can be used to parse a single region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7".
+             * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
+             * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
+             * If no precision is given, the one from the settings is used.
+             * 
+             */
+            static ParameterRegion parseRegion(
+                    std::string const& regionString,
+                    double precision=0.0);
+
+            /*
+             * Can be used to parse a vector of region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7;0.1<=p<=0.3,0.2<=q<=0.4".
+             * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
+             * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
+             * If no precision is given, the one from the settings is used.
+             * 
+             */
+            static std::vector<ParameterRegion> parseMultipleRegions(
+                    std::string const& regionsString,
+                    double precision=0.0);
+
+
+            /*
+             * Retrieves the regions that are specified in the settings.
+             * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
+             * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
+             * If no precision is given, the one from the settings is used.
+             * 
+             */
+            static std::vector<ParameterRegion> getRegionsFromSettings(double precision=0.0);
+
+            private:
+
+            std::map<VariableType, CoefficientType> const lowerBounds;
+            std::map<VariableType, CoefficientType> const upperBounds;
+            std::set<VariableType> variables;
+            RegionCheckResult checkResult;
+            std::map<VariableType, CoefficientType> satPoint;
+            std::map<VariableType, CoefficientType> violatedPoint;
+        };
+    }
+}
+
+#endif	/* STORM_MODELCHECKER_REGION_PARAMETERREGION_H */
+
diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp
index 7ed31bd37..3e2a46c7d 100644
--- a/src/modelchecker/region/SamplingModel.cpp
+++ b/src/modelchecker/region/SamplingModel.cpp
@@ -50,7 +50,7 @@ namespace storm {
             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<ParametricType, ConstantType>(parEntry.getValue())));
+                    samEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parEntry.getValue())));
                 }
                 else {
                     std::pair<ParametricType, ConstantType> searchedPair(parEntry.getValue(), storm::utility::zero<ConstantType>());
@@ -77,11 +77,7 @@ namespace storm {
             //write entries into evaluation table
             for(auto& tableEntry : this->evaluationTable){
                 tableEntry.second=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
-                        storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(
-                                tableEntry.first,
-                                point
-                            )
-                        );
+                        storm::utility::regions::evaluateFunction(tableEntry.first, point));
             }
             //write the instantiated values to the matrix according to the mapping
             for(auto& mappingPair : this->mapping){
diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h
index ffa50f221..0f5aabf22 100644
--- a/src/modelchecker/region/SamplingModel.h
+++ b/src/modelchecker/region/SamplingModel.h
@@ -8,7 +8,8 @@
 #ifndef STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H
 #define	STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H
 
-#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
+#include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
+#include "src/models/sparse/Dtmc.h"
 
 namespace storm {
     namespace modelchecker{
diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
new file mode 100644
index 000000000..e507a23d5
--- /dev/null
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
@@ -0,0 +1,734 @@
+#include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
+
+#include <chrono>
+
+#include "src/adapters/CarlAdapter.h"
+
+#include "src/modelchecker/region/ParameterRegion.h"
+#include "src/modelchecker/region/ApproximationModel.h"
+#include "src/modelchecker/region/SamplingModel.h"
+
+#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "src/utility/graph.h"
+#include "src/utility/vector.h"
+#include "src/utility/macros.h"
+
+#include "src/exceptions/InvalidPropertyException.h"
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/UnexpectedException.h"
+
+
+namespace storm {
+    namespace modelchecker {
+        
+        
+        template<typename ParametricType, typename ConstantType>
+        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc<ParametricType> const& model) : 
+                model(model),
+                eliminationModelChecker(model){
+            //intentionally left empty
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::canHandle(storm::logic::Formula const& formula) const {
+             //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ])
+            if(!formula.isStateFormula()){
+                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected a stateFormula");
+                return false;
+            }
+            if(!formula.asStateFormula().isProbabilityOperatorFormula()){
+                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected a probabilityOperatorFormula");
+                return false;
+            }
+            storm::logic::ProbabilityOperatorFormula const& probOpForm=formula.asStateFormula().asProbabilityOperatorFormula();
+            if(!probOpForm.hasBound()){
+                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: The formula has no bound");
+                return false;
+            }
+            if(!probOpForm.getSubformula().asPathFormula().isEventuallyFormula()){
+                STORM_LOG_DEBUG("Region Model Checker could not handle formula " << formula << " Reason: expected an eventually subformula");
+                return false;
+            }
+            if(model.getInitialStates().getNumberOfSetBits() != 1){
+                STORM_LOG_DEBUG("Input model is required to have exactly one initial state.");
+                return false;
+            }
+            return true;
+        }
+            
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::specifyFormula(storm::logic::Formula const& formula) {
+            std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now();
+            STORM_LOG_THROW(this->canHandle(formula), storm::exceptions::IllegalArgumentException, "Tried to specify a formula that can not be handled.");
+            
+            this->hasOnlyLinearFunctions=false;
+            this->isReachProbFunctionComputed=false;
+            this->isResultConstant=false;
+            this->smtSolver=nullptr;
+            
+            //Get subformula, target states
+            //Note: canHandle already ensures that the formula has the right shape and that the model has a single initial state.
+            this->probabilityOperatorFormula = std::unique_ptr<storm::logic::ProbabilityOperatorFormula>(new storm::logic::ProbabilityOperatorFormula(formula.asStateFormula().asProbabilityOperatorFormula()));
+            storm::logic::EventuallyFormula const& eventuallyFormula = this->probabilityOperatorFormula->getSubformula().asPathFormula().asEventuallyFormula();
+            std::unique_ptr<CheckResult> targetStatesResultPtr = this->eliminationModelChecker.check(eventuallyFormula.getSubformula());
+            storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector();
+            
+            // get a more simple model with a single target state, a single sink state and where states with constant outgoing transitions have been removed
+            // Note: also checks for linear functions and a constant result
+            computeSimplifiedModel(targetStates);
+            //now create the models used for sampling and approximation
+            this->samplingModel=std::make_shared<SamplingModel>(*this->simplifiedModel);
+            this->approximationModel=std::make_shared<ApproximationModel>(*this->simplifiedModel);
+           
+            //some information for statistics...
+            std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now();
+            this->timePreprocessing= timePreprocessingEnd - timePreprocessingStart;
+            this->numOfCheckedRegions=0;
+            this->numOfRegionsSolvedThroughSampling=0;
+            this->numOfRegionsSolvedThroughApproximation=0;
+            this->numOfRegionsSolvedThroughSubsystemSmt=0;
+            this->numOfRegionsSolvedThroughFullSmt=0;
+            this->numOfRegionsExistsBoth=0;
+            this->numOfRegionsAllSat=0;
+            this->numOfRegionsAllViolated=0;
+            this->timeCheckRegion=std::chrono::high_resolution_clock::duration::zero();
+            this->timeSampling=std::chrono::high_resolution_clock::duration::zero();
+            this->timeApproximation=std::chrono::high_resolution_clock::duration::zero();
+            this->timeMDPBuild=std::chrono::high_resolution_clock::duration::zero();
+            this->timeSubsystemSmt=std::chrono::high_resolution_clock::duration::zero();
+            this->timeFullSmt=std::chrono::high_resolution_clock::duration::zero();
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::computeSimplifiedModel(storm::storage::BitVector const& targetStates){
+                        
+            //Compute the subset of states that have a probability of 0 or 1, respectively and reduce the considered states accordingly.
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->model, storm::storage::BitVector(this->model.getNumberOfStates(),true), targetStates);
+            storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
+            storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
+            storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+            // If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction.
+            if (this->model.getInitialStates().isDisjointFrom(maybeStates)) {
+                STORM_LOG_WARN("The probability of the initial state is constant (0 or 1)");
+                this->reachProbFunction = statesWithProbability0.get(*(this->model.getInitialStates().begin())) ? storm::utility::zero<ParametricType>() : storm::utility::one<ParametricType>();
+                this->isReachProbFunctionComputed=true;
+                this->isResultConstant=true;
+                this->hasOnlyLinearFunctions=true;
+            }
+            // Determine the set of states that is reachable from the initial state without jumping over a target state.
+            storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->model.getTransitionMatrix(), this->model.getInitialStates(), maybeStates, statesWithProbability1);
+            // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
+            maybeStates &= reachableStates;
+            // Create a vector for the probabilities to go to a state with probability 1 in one step.
+            std::vector<ParametricType> oneStepProbabilities = this->model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
+            // Determine the initial state of the sub-model.
+            storm::storage::sparse::state_type initialState = *(this->model.getInitialStates() % maybeStates).begin();
+            // We then build the submatrix that only has the transitions of the maybe states.
+            storm::storage::SparseMatrix<ParametricType> submatrix = this->model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
+            
+            // eliminate all states with only constant outgoing transitions
+            //TODO: maybe also states with constant incoming tranistions. THEN the ordering of the eliminated states does matter.
+            std::chrono::high_resolution_clock::time_point timeInitialStateEliminationStart = std::chrono::high_resolution_clock::now();
+            // Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
+            auto flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix);
+            auto flexibleBackwardTransitions= this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix.transpose(), true);
+            // Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated.
+            storm::storage::BitVector subsystem = storm::storage::BitVector(submatrix.getRowCount(), true);
+            //temporarily unselect the initial state to skip it.
+            subsystem.set(initialState, false);
+            this->hasOnlyLinearFunctions=true;
+            bool allReachableFunctionsAreConstant=true;
+            boost::optional<std::vector<ParametricType>> missingStateRewards;
+            for (auto const& state : subsystem) {
+                bool stateHasOnlyConstantOutgoingTransitions=true;
+                for(auto const& entry : submatrix.getRow(state)){
+                    if(!this->parametricTypeComparator.isConstant(entry.getValue())){
+                        allReachableFunctionsAreConstant=false;
+                        stateHasOnlyConstantOutgoingTransitions=false;
+                        this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(entry.getValue());
+                    }
+                }
+                if(!this->parametricTypeComparator.isConstant(oneStepProbabilities[state])){
+                    allReachableFunctionsAreConstant=false;
+                    stateHasOnlyConstantOutgoingTransitions=false;
+                    this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(oneStepProbabilities[state]);
+                }
+                if(stateHasOnlyConstantOutgoingTransitions){
+                    this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
+                    subsystem.set(state,false);
+                }
+            }
+            subsystem.set(initialState, true);
+            std::chrono::high_resolution_clock::time_point timeInitialStateEliminationEnd = std::chrono::high_resolution_clock::now();
+            this->timeInitialStateElimination = timeInitialStateEliminationEnd-timeInitialStateEliminationStart;
+            STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl);
+            std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl;
+            
+            if(allReachableFunctionsAreConstant){
+                //Check if this is also the case for the initial state
+                for(auto const& entry : submatrix.getRow(initialState)){
+                    allReachableFunctionsAreConstant&=this->parametricTypeComparator.isConstant(entry.getValue());
+                }
+                allReachableFunctionsAreConstant&=this->parametricTypeComparator.isConstant(oneStepProbabilities[initialState]);
+                // Set the flag accordingly
+                this->isResultConstant=allReachableFunctionsAreConstant;
+                STORM_LOG_WARN_COND(!this->isResultConstant, "For the given property, the reachability probability is constant, i.e., independent of the region");
+            }
+            //Build the simplified model
+            //The matrix. The flexibleTransitions matrix might have empty rows where states have been eliminated.
+            //The new matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
+            std::vector<storm::storage::sparse::state_type> newStateIndexMap(flexibleTransitions.getNumberOfRows(), flexibleTransitions.getNumberOfRows()); //initialize with some illegal index to easily check if a transition leads to an unselected state
+            storm::storage::sparse::state_type newStateIndex=0;
+            for(auto const& state : subsystem){
+                newStateIndexMap[state]=newStateIndex;
+                ++newStateIndex;
+            }
+            //We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead
+            storm::storage::sparse::state_type numStates=newStateIndex+2;
+            storm::storage::sparse::state_type targetState=numStates-2;
+            storm::storage::sparse::state_type sinkState= numStates-1;
+            //We can now fill in the data.
+            storm::storage::SparseMatrixBuilder<ParametricType> matrixBuilder(numStates, numStates);
+            for(storm::storage::sparse::state_type oldStateIndex : subsystem){ 
+                ParametricType valueToSinkState=storm::utility::regions::getNewFunction<ParametricType, CoefficientType>(storm::utility::one<CoefficientType>());
+                //go through columns:
+                for(auto const& entry: flexibleTransitions.getRow(oldStateIndex)){ 
+                    STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=flexibleTransitions.getNumberOfRows(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated.");
+                    valueToSinkState-=entry.getValue();
+                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],entry.getValue());
+                }
+                //transition to target state
+                if(!this->parametricTypeComparator.isZero(oneStepProbabilities[oldStateIndex])){ 
+                    valueToSinkState-=oneStepProbabilities[oldStateIndex];
+                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, oneStepProbabilities[oldStateIndex]);
+                }
+                //transition to sink state
+                if(!this->parametricTypeComparator.isZero(valueToSinkState)){ 
+                    matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, valueToSinkState);
+                }
+            }
+            //add self loops on the additional states (i.e., target and sink)
+            matrixBuilder.addNextValue(targetState, targetState, storm::utility::one<ParametricType>());
+            matrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one<ParametricType>());
+            //The labeling
+            storm::models::sparse::StateLabeling labeling(numStates);
+            storm::storage::BitVector initLabel(numStates, false);
+            initLabel.set(newStateIndexMap[initialState], true);
+            labeling.addLabel("init", std::move(initLabel));
+            storm::storage::BitVector targetLabel(numStates, false);
+            targetLabel.set(targetState, true);
+            labeling.addLabel("target", std::move(targetLabel));
+            storm::storage::BitVector sinkLabel(numStates, false);
+            sinkLabel.set(sinkState, true);
+            labeling.addLabel("sink", std::move(sinkLabel));
+            // other ingredients
+            boost::optional<std::vector<ParametricType>> noStateRewards;
+            boost::optional<storm::storage::SparseMatrix<ParametricType>> noTransitionRewards;  
+            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;            
+            // the final model
+            this->simplifiedModel = std::make_shared<storm::models::sparse::Dtmc<ParametricType>>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
+        }
+                    
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkRegion(ParameterRegion& region) {
+            std::chrono::high_resolution_clock::time_point timeCheckRegionStart = std::chrono::high_resolution_clock::now();
+            ++this->numOfCheckedRegions;
+            
+            STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to analyze a region although no property has been specified" );
+            STORM_LOG_DEBUG("Analyzing the region " << region.toString());
+            
+            
+            //switches for the different steps.
+            bool done=false;
+            bool doApproximation=this->hasOnlyLinearFunctions; // this approach is only correct if the model has only linear functions
+            bool doSampling=true;
+            bool doSubsystemSmt=false; //this->hasOnlyLinearFunctions;  // this approach is only correct if the model has only linear functions
+            bool doFullSmt=true;
+            
+            if(!done && this->isResultConstant){
+                STORM_LOG_DEBUG("Checking a region although the result is constant, i.e., independent of the region. This makes sense none.");
+                STORM_LOG_THROW(this->parametricTypeComparator.isConstant(getReachProbFunction()), storm::exceptions::UnexpectedException, "The result was assumed to be constant but it isn't.");
+                if(valueIsInBoundOfFormula(storm::utility::regions::getConstantPart(getReachProbFunction()))){
+                    region.setCheckResult(RegionCheckResult::ALLSAT);
+                }
+                else{
+                    region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                }
+                done=true;
+            }
+            
+            std::chrono::high_resolution_clock::time_point timeApproximationStart = std::chrono::high_resolution_clock::now();
+            std::vector<ConstantType> lowerBounds;
+            std::vector<ConstantType> upperBounds;
+            if(!done && doApproximation){
+                STORM_LOG_DEBUG("Checking approximative probabilities...");
+                if (region.getCheckResult()==RegionCheckResult::UNKNOWN){
+                    //Sample a single point to know whether we should try to prove ALLSAT or ALLVIOLATED
+                    checkPoint(region,region.getLowerBounds(), false);
+                }
+                if(checkApproximativeProbabilities(region, lowerBounds, upperBounds)){
+                    ++this->numOfRegionsSolvedThroughApproximation;
+                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through approximation.");
+                    done=true;
+                }
+            }
+            std::chrono::high_resolution_clock::time_point timeApproximationEnd = std::chrono::high_resolution_clock::now();
+            
+            std::chrono::high_resolution_clock::time_point timeSamplingStart = std::chrono::high_resolution_clock::now();
+            if(!done && doSampling){
+                STORM_LOG_DEBUG("Checking sample points...");
+                if(checkSamplePoints(region)){
+                    ++this->numOfRegionsSolvedThroughSampling;
+                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through sampling.");
+                    done=true;
+                }
+            }
+            std::chrono::high_resolution_clock::time_point timeSamplingEnd = std::chrono::high_resolution_clock::now();
+            
+            std::chrono::high_resolution_clock::time_point timeSubsystemSmtStart = std::chrono::high_resolution_clock::now();
+            if(!done && doSubsystemSmt){
+                STORM_LOG_WARN("SubsystemSmt approach not yet implemented");
+            }
+            std::chrono::high_resolution_clock::time_point timeSubsystemSmtEnd = std::chrono::high_resolution_clock::now();
+            
+            std::chrono::high_resolution_clock::time_point timeFullSmtStart = std::chrono::high_resolution_clock::now();
+            if(!done && doFullSmt){
+                STORM_LOG_DEBUG("Checking with Smt Solving...");
+                if(checkFullSmt(region)){
+                    ++this->numOfRegionsSolvedThroughFullSmt;
+                    STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through Smt Solving.");
+                    done=true;
+                }
+            }
+            std::chrono::high_resolution_clock::time_point timeFullSmtEnd = std::chrono::high_resolution_clock::now();
+            
+            //some information for statistics...
+            std::chrono::high_resolution_clock::time_point timeCheckRegionEnd = std::chrono::high_resolution_clock::now();
+            this->timeCheckRegion += timeCheckRegionEnd-timeCheckRegionStart;
+            this->timeSampling += timeSamplingEnd - timeSamplingStart;
+            this->timeApproximation += timeApproximationEnd - timeApproximationStart;
+            this->timeSubsystemSmt += timeSubsystemSmtEnd - timeSubsystemSmtStart;
+            this->timeFullSmt += timeFullSmtEnd - timeFullSmtStart;
+            switch(region.getCheckResult()){
+                case RegionCheckResult::EXISTSBOTH:
+                    ++this->numOfRegionsExistsBoth;
+                    break;
+                case RegionCheckResult::ALLSAT:
+                    ++this->numOfRegionsAllSat;
+                    break;
+                case RegionCheckResult::ALLVIOLATED:
+                    ++this->numOfRegionsAllViolated;
+                    break;
+                default:
+                    break;
+            }
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkRegions(std::vector<ParameterRegion>& regions) {
+            STORM_LOG_DEBUG("Checking " << regions.size() << "regions.");
+            std::cout << "Checking " << regions.size() << " regions. Progress: ";
+            std::cout.flush();
+                    
+            uint_fast64_t progress=0;
+            uint_fast64_t checkedRegions=0;
+            for(auto& region : regions){
+                this->checkRegion(region);
+                if((checkedRegions++)*10/regions.size()==progress){
+                    std::cout << progress++;
+                    std::cout.flush();
+                }
+            }
+            std::cout << " done!" << std::endl;
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkApproximativeProbabilities(ParameterRegion& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds) {
+            STORM_LOG_THROW(this->hasOnlyLinearFunctions, storm::exceptions::UnexpectedException, "Tried to perform approximative method (only applicable if all functions are linear), but there are nonlinear functions.");
+            std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now();
+            this->approximationModel->instantiate(region);
+            std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now();
+            this->timeMDPBuild += timeMDPBuildEnd-timeMDPBuildStart;
+            
+            // Decide whether we should minimize or maximize and whether to prove allsat or allviolated. (Hence, there are 4 cases)
+            bool formulaHasUpperBound = this->probabilityOperatorFormula->getComparisonType()==storm::logic::ComparisonType::Less || this->probabilityOperatorFormula->getComparisonType()==storm::logic::ComparisonType::LessEqual;
+            STORM_LOG_THROW((formulaHasUpperBound != (this->probabilityOperatorFormula->getComparisonType()==storm::logic::ComparisonType::Greater || this->probabilityOperatorFormula->getComparisonType()==storm::logic::ComparisonType::GreaterEqual)),
+                    storm::exceptions::UnexpectedException, "Unexpected comparison Type of formula");
+            switch (region.getCheckResult()){
+                case RegionCheckResult::EXISTSSAT:
+                case RegionCheckResult::UNKNOWN: 
+                    //Try to prove ALLSAT
+                    if(formulaHasUpperBound){ 
+                        upperBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Maximize);
+                        lowerBounds=std::vector<ConstantType>();
+                        if(valueIsInBoundOfFormula(upperBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                            region.setCheckResult(RegionCheckResult::ALLSAT);
+                            return true;
+                        }
+                    }
+                    else{
+                        lowerBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Minimize);
+                        upperBounds=std::vector<ConstantType>();
+                        if(valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                            region.setCheckResult(RegionCheckResult::ALLSAT);
+                            return true;
+                        }
+                    }
+                    break;
+                case RegionCheckResult::EXISTSVIOLATED:
+                    //Try to prove ALLVIOLATED
+                    if(formulaHasUpperBound){ 
+                        lowerBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Minimize);
+                        upperBounds=std::vector<ConstantType>();
+                        if(!valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                            return true;
+                        }
+                    }
+                    else{
+                        upperBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Maximize);
+                        lowerBounds=std::vector<ConstantType>();
+                        if(!valueIsInBoundOfFormula(upperBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                            return true;
+                        }
+                    }
+                    break;
+                default:
+                    STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The checkresult of the current region should not be conclusive, i.e. it should be either EXISTSSAT or EXISTSVIOLATED or UNKNOWN");
+            }
+            
+            if(region.getCheckResult()==RegionCheckResult::UNKNOWN){
+                //In this case, it makes sense to try the other optimality Type. Again, 4 cases (Although the ALLSAT cases should not be reachable, since we already tried to prove this above).
+                if(lowerBounds.empty()){
+                    lowerBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Minimize);
+                    if(!formulaHasUpperBound && valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                        region.setCheckResult(RegionCheckResult::ALLSAT);
+                        return true;
+                    }
+                    if(formulaHasUpperBound && !valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                        region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                        return true;
+                    }
+                }
+                if(upperBounds.empty()){
+                    upperBounds = this->approximationModel->computeReachabilityProbabilities(storm::logic::OptimalityType::Maximize);
+                    if(formulaHasUpperBound && valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                        region.setCheckResult(RegionCheckResult::ALLSAT);
+                        return true;
+                    }
+                    if(!formulaHasUpperBound && !valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
+                        region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                        return true;
+                    }
+                }
+            }
+            //if we reach this point than the result is still inconclusive.
+            return false;            
+        }
+        
+              template<typename ParametricType, typename ConstantType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkSamplePoints(ParameterRegion& region) {
+            auto samplingPoints = region.getVerticesOfRegion(region.getVariables()); //test the 4 corner points
+            for (auto const& point : samplingPoints){
+                if(checkPoint(region, point, false)){
+                    return true;
+                }            
+            }
+            return false;
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkPoint(ParameterRegion& region, std::map<VariableType, CoefficientType>const& point, bool viaReachProbFunction) {
+            bool valueInBoundOfFormula;
+            if(viaReachProbFunction){
+                valueInBoundOfFormula = this->valueIsInBoundOfFormula(storm::utility::regions::evaluateFunction(getReachProbFunction(), point));
+            }
+            else{
+                this->samplingModel->instantiate(point);
+                valueInBoundOfFormula=this->valueIsInBoundOfFormula(this->samplingModel->computeReachabilityProbabilities()[*this->samplingModel->getModel()->getInitialStates().begin()]);
+            }
+                
+            if(valueInBoundOfFormula){
+                if (region.getCheckResult()!=RegionCheckResult::EXISTSSAT){
+                    region.setSatPoint(point);
+                    if(region.getCheckResult()==RegionCheckResult::EXISTSVIOLATED){
+                        region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                        return true;
+                    }
+                    region.setCheckResult(RegionCheckResult::EXISTSSAT);
+                }
+            }
+            else{
+                if (region.getCheckResult()!=RegionCheckResult::EXISTSVIOLATED){
+                    region.setViolatedPoint(point);
+                    if(region.getCheckResult()==RegionCheckResult::EXISTSSAT){
+                        region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                        return true;
+                    }
+                    region.setCheckResult(RegionCheckResult::EXISTSVIOLATED);
+                }
+            }
+            return false;
+        }
+         
+        template<typename ParametricType, typename ConstantType>
+        ParametricType SparseDtmcRegionModelChecker<ParametricType, ConstantType>::getReachProbFunction() {
+            if(!this->isReachProbFunctionComputed){
+                std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionStart = std::chrono::high_resolution_clock::now();
+                //get the one step probabilities and the transition matrix of the simplified model without target/sink state
+                //TODO check if this takes long... we could also store the oneSteps while creating the simplified model. Or(?) we keep the matrix the way it is and give the target state one step probability 1.
+                storm::storage::SparseMatrix<ParametricType> backwardTransitions(this->simplifiedModel->getBackwardTransitions());
+                std::vector<ParametricType> oneStepProbabilities(this->simplifiedModel->getNumberOfStates()-2, storm::utility::zero<ParametricType>());
+                for(auto const& entry : backwardTransitions.getRow(*(this->simplifiedModel->getStates("target").begin()))){
+                    if(entry.getColumn()<oneStepProbabilities.size()){
+                        oneStepProbabilities[entry.getColumn()]=entry.getValue();
+                    } //else case: only holds for the entry that corresponds to the selfloop on the target state..
+                }
+                storm::storage::BitVector maybeStates=~(this->simplifiedModel->getStates("target") | this->simplifiedModel->getStates("sink"));
+                backwardTransitions=backwardTransitions.getSubmatrix(false,maybeStates,maybeStates);
+                storm::storage::SparseMatrix<ParametricType> forwardTransitions=this->simplifiedModel->getTransitionMatrix().getSubmatrix(false,maybeStates,maybeStates);
+                //now compute the functions using methods from elimination model checker
+                storm::storage::BitVector newInitialStates = this->simplifiedModel->getInitialStates() % maybeStates;
+                storm::storage::BitVector phiStates(this->simplifiedModel->getNumberOfStates(), true);
+                boost::optional<std::vector<ParametricType>> noStateRewards;
+                std::vector<std::size_t> statePriorities = this->eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities);
+                this->reachProbFunction=this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, this->simplifiedModel->getStates("target"), noStateRewards, statePriorities);
+                std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionEnd = std::chrono::high_resolution_clock::now();
+            /*    std::string funcStr = " (/ " +
+                                    this->reachProbFunction.nominator().toString(false, true) + " " +
+                                    this->reachProbFunction.denominator().toString(false, true) +
+                                " )";
+                std::cout << std::endl <<"the resulting reach prob function is " << std::endl << funcStr << std::endl << std::endl;*/
+                STORM_LOG_DEBUG("Computed reachProbFunction");
+                this->isReachProbFunctionComputed=true;
+                this->timeComputeReachProbFunction = timeComputeReachProbFunctionEnd-timeComputeReachProbFunctionStart;
+            }
+            return this->reachProbFunction;
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkFullSmt(ParameterRegion& region) {
+            if (region.getCheckResult()==RegionCheckResult::UNKNOWN){
+                //Sampling needs to be done (on a single point)
+                checkPoint(region,region.getLowerBounds(), true);
+            }
+            
+            if(this->smtSolver==nullptr){
+                initializeSMTSolver(this->smtSolver, getReachProbFunction(),*this->probabilityOperatorFormula);
+            }
+            
+            this->smtSolver->push();
+            
+            //add constraints for the region
+            for(auto const& variable : region.getVariables()) {
+                storm::utility::regions::addParameterBoundsToSmtSolver(this->smtSolver, variable, storm::logic::ComparisonType::GreaterEqual, region.getLowerBound(variable));
+                storm::utility::regions::addParameterBoundsToSmtSolver(this->smtSolver, variable, storm::logic::ComparisonType::LessEqual, region.getUpperBound(variable));
+            }
+            
+            //add constraint that states what we want to prove            
+            VariableType proveAllSatVar=storm::utility::regions::getVariableFromString<VariableType>("storm_proveAllSat");     
+            VariableType proveAllViolatedVar=storm::utility::regions::getVariableFromString<VariableType>("storm_proveAllViolated");     
+            switch(region.getCheckResult()){
+                case RegionCheckResult::EXISTSBOTH:
+                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::EXISTSBOTH), "checkFullSmt invoked although the result is already clear (EXISTSBOTH). Will validate this now...");
+                case RegionCheckResult::ALLSAT:
+                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::ALLSAT), "checkFullSmt invoked although the result is already clear (ALLSAT). Will validate this now...");
+                case RegionCheckResult::EXISTSSAT:
+                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllSatVar, true);
+                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllViolatedVar, false);
+                    break;
+                case RegionCheckResult::ALLVIOLATED:
+                    STORM_LOG_WARN_COND((region.getCheckResult()!=RegionCheckResult::ALLVIOLATED), "checkFullSmt invoked although the result is already clear (ALLVIOLATED). Will validate this now...");
+                case RegionCheckResult::EXISTSVIOLATED:
+                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllSatVar, false);
+                    storm::utility::regions::addBoolVariableToSmtSolver(this->smtSolver, proveAllViolatedVar, true);
+                    break;
+                default:
+                STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not handle the current region CheckResult: " << region.checkResultToString());
+            }
+            
+            storm::solver::SmtSolver::CheckResult solverResult= this->smtSolver->check();
+            this->smtSolver->pop();
+            
+            switch(solverResult){
+                case storm::solver::SmtSolver::CheckResult::Sat:
+                    switch(region.getCheckResult()){
+                        case RegionCheckResult::EXISTSSAT:
+                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                            //There is also a violated point
+                            STORM_LOG_WARN("Extracting a violated point from the smt solver is not yet implemented!");
+                            break;
+                        case RegionCheckResult::EXISTSVIOLATED:
+                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                            //There is also a sat point
+                            STORM_LOG_WARN("Extracting a sat point from the smt solver is not yet implemented!");
+                            break;
+                        case RegionCheckResult::EXISTSBOTH:
+                            //That was expected
+                            STORM_LOG_WARN("result EXISTSBOTH Validated!");
+                            break;
+                        default:
+                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The solver gave an unexpected result (sat)");
+                    }
+                    return true;
+                case storm::solver::SmtSolver::CheckResult::Unsat:
+                    switch(region.getCheckResult()){
+                        case RegionCheckResult::EXISTSSAT:
+                            region.setCheckResult(RegionCheckResult::ALLSAT);
+                            break;
+                        case RegionCheckResult::EXISTSVIOLATED:
+                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                            break;
+                        case RegionCheckResult::ALLSAT:
+                            //That was expected...
+                            STORM_LOG_WARN("result ALLSAT Validated!");
+                            break;
+                        case RegionCheckResult::ALLVIOLATED:
+                            //That was expected...
+                            STORM_LOG_WARN("result ALLVIOLATED Validated!");
+                            break;
+                        default:
+                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The solver gave an unexpected result (unsat)");
+                    }
+                    return true;
+                case storm::solver::SmtSolver::CheckResult::Unknown:
+                default:
+                    STORM_LOG_WARN("The SMT solver was not able to compute a result for this region. (Timeout? Memout?)");
+                    if(this->smtSolver->isNeedsRestart()){
+                        initializeSMTSolver(this->smtSolver,getReachProbFunction(), *this->probabilityOperatorFormula);
+                    }
+                    return false;
+            }
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, const ParametricType& reachProbFunc, const storm::logic::ProbabilityOperatorFormula& formula) {
+                    
+            storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions..
+            solver = std::shared_ptr<storm::solver::Smt2SmtSolver>(new storm::solver::Smt2SmtSolver(manager, true));
+            
+            ParametricType bound= storm::utility::regions::convertNumber<double, ParametricType>(this->probabilityOperatorFormula->getBound());
+            
+            // To prove that the property is satisfied in the initial state for all parameters,
+            // we ask the solver whether the negation of the property is satisfiable and invert the answer.
+            // In this case, assert that this variable is true:
+            VariableType proveAllSatVar=storm::utility::regions::getNewVariable<VariableType>("storm_proveAllSat", storm::utility::regions::VariableSort::VS_BOOL);
+            
+            //Example:
+            //Property:    P<=p [ F 'target' ] holds iff...
+            // f(x)         <= p
+            // Hence: If  f(x) > p is unsat, the property is satisfied for all parameters.
+            
+            storm::logic::ComparisonType proveAllSatRel; //the relation from the property needs to be inverted
+            switch (this->probabilityOperatorFormula->getComparisonType()) {
+                case storm::logic::ComparisonType::Greater:
+                    proveAllSatRel=storm::logic::ComparisonType::LessEqual;
+                    break;
+                case storm::logic::ComparisonType::GreaterEqual:
+                    proveAllSatRel=storm::logic::ComparisonType::Less;
+                    break;
+                case storm::logic::ComparisonType::Less:
+                    proveAllSatRel=storm::logic::ComparisonType::GreaterEqual;
+                    break;
+                case storm::logic::ComparisonType::LessEqual:
+                    proveAllSatRel=storm::logic::ComparisonType::Greater;
+                    break;
+                default:
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
+            }
+            storm::utility::regions::addGuardedConstraintToSmtSolver(solver, proveAllSatVar, reachProbFunc, proveAllSatRel, bound);
+            
+            // To prove that the property is violated in the initial state for all parameters,
+            // we ask the solver whether the the property is satisfiable and invert the answer.
+            // In this case, assert that this variable is true:
+            VariableType proveAllViolatedVar=storm::utility::regions::getNewVariable<VariableType>("storm_proveAllViolated", storm::utility::regions::VariableSort::VS_BOOL);        
+            
+            //Example:
+            //Property:    P<=p [ F 'target' ] holds iff...
+            // f(x)         <= p
+            // Hence: If f(x)  <= p is unsat, the property is violated for all parameters. 
+            storm::logic::ComparisonType proveAllViolatedRel = this->probabilityOperatorFormula->getComparisonType();
+            storm::utility::regions::addGuardedConstraintToSmtSolver(solver, proveAllViolatedVar, reachProbFunc, proveAllViolatedRel, bound);          
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        template<typename ValueType>
+        bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::valueIsInBoundOfFormula(ValueType value){
+            STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to compare a value to the bound of a formula, but no formula specified.");
+            double valueAsDouble = storm::utility::regions::convertNumber<ValueType, double>(value);
+            switch (this->probabilityOperatorFormula->getComparisonType()) {
+                case storm::logic::ComparisonType::Greater:
+                    return (valueAsDouble > this->probabilityOperatorFormula->getBound());
+                case storm::logic::ComparisonType::GreaterEqual:
+                    return (valueAsDouble >= this->probabilityOperatorFormula->getBound());
+                case storm::logic::ComparisonType::Less:
+                    return (valueAsDouble < this->probabilityOperatorFormula->getBound());
+                case storm::logic::ComparisonType::LessEqual:
+                    return (valueAsDouble <= this->probabilityOperatorFormula->getBound());
+                default:
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
+            }
+        }
+
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::printStatisticsToStream(std::ostream& outstream) {
+            
+            if(this->probabilityOperatorFormula==nullptr){
+                outstream << "Region Model Checker Statistics Error: No formula specified." << std::endl; 
+                return;
+            }
+            
+            std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timePreprocessing);
+            std::chrono::milliseconds timeInitialStateEliminationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeInitialStateElimination);
+            std::chrono::milliseconds timeComputeReachProbFunctionInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeComputeReachProbFunction);
+            std::chrono::milliseconds timeCheckRegionInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeCheckRegion);
+            std::chrono::milliseconds timeSammplingInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeSampling);
+            std::chrono::milliseconds timeApproximationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeApproximation);
+            std::chrono::milliseconds timeMDPBuildInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeMDPBuild);
+            std::chrono::milliseconds timeFullSmtInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(this->timeFullSmt);
+            
+            std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ...
+            std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeOverall);
+            
+            uint_fast64_t numOfSolvedRegions= this->numOfRegionsExistsBoth + this->numOfRegionsAllSat + this->numOfRegionsAllViolated;
+            
+            outstream << std::endl << "Region Model Checker Statistics:" << std::endl;
+            outstream << "-----------------------------------------------" << std::endl;
+            outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl;
+            outstream << "Simplified model: " << this->simplifiedModel->getNumberOfStates() << " states, " << this->simplifiedModel->getNumberOfTransitions() << " transitions" << std::endl;
+            outstream << "Formula: " << *this->probabilityOperatorFormula << std::endl;
+            outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl;
+            outstream << "Number of checked regions: " << this->numOfCheckedRegions << std::endl;
+            outstream << "  Number of solved regions:  " <<  numOfSolvedRegions << "(" << numOfSolvedRegions*100/this->numOfCheckedRegions << "%)" <<  std::endl;
+            outstream << "    AllSat:      " <<  this->numOfRegionsAllSat << "(" << this->numOfRegionsAllSat*100/this->numOfCheckedRegions << "%)" <<  std::endl;
+            outstream << "    AllViolated: " <<  this->numOfRegionsAllViolated << "(" << this->numOfRegionsAllViolated*100/this->numOfCheckedRegions << "%)" <<  std::endl;
+            outstream << "    ExistsBoth:  " <<  this->numOfRegionsExistsBoth << "(" << this->numOfRegionsExistsBoth*100/this->numOfCheckedRegions << "%)" <<  std::endl;
+            outstream << "    Unsolved:    " <<  this->numOfCheckedRegions - numOfSolvedRegions << "(" << (this->numOfCheckedRegions - numOfSolvedRegions)*100/this->numOfCheckedRegions << "%)" <<  std::endl;
+            outstream << "  --  Note: %-numbers are relative to the NUMBER of regions, not the size of their area --" <<  std::endl;
+            outstream << "  " << this->numOfRegionsSolvedThroughSampling << " regions solved through Sampling" << std::endl;
+            outstream << "  " << this->numOfRegionsSolvedThroughApproximation << " regions solved through Approximation" << std::endl;
+            outstream << "  " << this->numOfRegionsSolvedThroughSubsystemSmt << " regions solved through SubsystemSmt" << std::endl;
+            outstream << "  " << this->numOfRegionsSolvedThroughFullSmt << " regions solved through FullSmt" << std::endl;
+            outstream << std::endl;
+            outstream << "Running times:" << std::endl;
+            outstream << "  " << timeOverallInMilliseconds.count() << "ms overall (excluding model parsing)" << std::endl;
+            outstream << "  " << timePreprocessingInMilliseconds.count() << "ms Preprocessing including... " << std::endl;
+            outstream << "    " << timeInitialStateEliminationInMilliseconds.count() << "ms Initial state elimination of const transitions" << std::endl;
+            outstream << "  " << timeComputeReachProbFunctionInMilliseconds.count() << "ms to compute the reachability probability function" << std::endl;
+            outstream << "  " << timeCheckRegionInMilliseconds.count() << "ms Region Check including... " << std::endl;
+            outstream << "    " << timeSammplingInMilliseconds.count() << "ms Sampling " << std::endl;
+            outstream << "    " << timeApproximationInMilliseconds.count() << "ms Approximation including... " << std::endl;
+            outstream << "      " << timeMDPBuildInMilliseconds.count() << "ms to build the MDP" << std::endl;
+            outstream << "    " << timeFullSmtInMilliseconds.count() << "ms Full Smt solving" << std::endl;
+            outstream << "-----------------------------------------------" << std::endl;
+            
+        }
+        
+#ifdef STORM_HAVE_CARL
+        template class SparseDtmcRegionModelChecker<storm::RationalFunction, double>;
+#endif
+        //note: for other template instantiations, add rules for the typedefs of VariableType and CoefficientType in utility/regions.h
+        
+    } // namespace modelchecker
+} // namespace storm
diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.h b/src/modelchecker/region/SparseDtmcRegionModelChecker.h
new file mode 100644
index 000000000..0dca1f257
--- /dev/null
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.h
@@ -0,0 +1,205 @@
+#ifndef STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_
+#define STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_
+
+#include "src/storage/sparse/StateType.h"
+#include "src/models/sparse/Dtmc.h"
+#include "src/utility/constants.h"
+#include "src/utility/regions.h"
+#include "src/solver/Smt2SmtSolver.h"
+#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
+
+namespace storm {
+    namespace modelchecker {
+        
+        template<typename ParametricType, typename ConstantType>
+        class SparseDtmcRegionModelChecker {
+        public:
+            
+            typedef typename storm::utility::regions::VariableType<ParametricType> VariableType;
+            typedef typename storm::utility::regions::CoefficientType<ParametricType> CoefficientType;
+            
+            /*!
+             * The possible results for a single Parameter region
+             */
+            enum class RegionCheckResult { 
+                UNKNOWN, /*!< the result is unknown */
+                EXISTSSAT, /*!< the formula is satisfied for at least one parameter evaluation that lies in the given region */
+                EXISTSVIOLATED, /*!< the formula is violated for at least one parameter evaluation that lies in the given region */
+                EXISTSBOTH, /*!< the formula is satisfied for some parameters but also violated for others */
+                ALLSAT, /*!< the formula is satisfied for all parameters in the given region */
+                ALLVIOLATED /*!< the formula is violated for all parameters in the given region */
+            };
+            
+            class ParameterRegion;
+            class ApproximationModel;
+            class SamplingModel;
+            
+            explicit SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc<ParametricType> const& model);
+
+            /*!
+             * Checks if the given formula can be handled by This region model checker
+             * @param formula the formula to be checked
+             */
+            bool canHandle(storm::logic::Formula const& formula) const;
+            
+            /*!
+             * Specifies the considered formula.
+             * A few preprocessing steps are performed.
+             * If another formula has been specified before, all 'context' regarding the old formula is lost.
+             * 
+             * @param formula the formula to be considered.
+             */
+            void specifyFormula(storm::logic::Formula const& formula);
+
+            /*!
+             * Checks whether the given formula holds for all parameters that lie in the given region.
+             * Sets the region checkresult accordingly. Moreover, region.satPoint and/or an region.violatedPoint will be set.
+             * 
+             * @note A formula has to be specified first.
+             * 
+             * @param region The considered region
+             * 
+             */
+            void checkRegion(ParameterRegion& region);
+            
+            /*!
+             * Checks for every given region whether the specified formula holds for all parameters that lie in that region.
+             * Sets the region checkresult accordingly. Moreover, region.satPoint and/or an region.violatedPoint will be set.
+             * 
+             * @note A formula has to be specified first.
+             * 
+             * @param region The considered region
+             */
+            void checkRegions(std::vector<ParameterRegion>& regions);
+            
+            /*!
+             * Prints statistical information to the given stream.
+             */
+            void printStatisticsToStream(std::ostream& outstream);
+            
+            
+            
+        private:
+            
+            /*!
+             * Computes a model with a single target and a single sink state.
+             * Eliminates all states for which the outgoing transitions are constant.
+             * 
+             * Also checks whether the non constant functions are linear and whether the model checking result is independent of parameters (i.e., constant)
+             * The flags of This model checker  are set accordingly.
+             */
+            void computeSimplifiedModel(storm::storage::BitVector const& targetStates);
+            
+            /*!
+             * Instantiates the approximation model to compute bounds on the maximal/minimal reachability probability.
+             * If the current region result is EXISTSSAT (or EXISTSVIOLATED), then this function tries to prove ALLSAT (or ALLVIOLATED).
+             * If this succeeded, then the region check result is changed accordingly.
+             * If the current region result is UNKNOWN, then this function first tries to prove ALLSAT and if that failed, it tries to prove ALLVIOLATED.
+             * In any case, the computed bounds are written to the given lowerBounds/upperBounds.
+             * However, if only the lowerBounds (or upperBounds) have been computed, the other vector is set to a vector of size 0.
+             * True is returned iff either ALLSAT or ALLVIOLATED could be proved.
+             */
+            bool checkApproximativeProbabilities(ParameterRegion& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds); 
+            
+            /*!
+             * Checks the value of the function at some sampling points within the given region.
+             * May set the satPoint and violatedPoint of the regions if they are not yet specified and such points are found
+             * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
+             * 
+             * @return true if an violated point as well as a sat point has been found during the process
+             */
+            bool checkSamplePoints(ParameterRegion& region);
+            
+            /*!
+             * Checks the value of the function at the given sampling point.
+             * May set the satPoint and violatedPoint of the regions if thy are not yet specified and such point is given.
+             * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
+             * 
+             * @param viaReachProbFunction if set, the sampling will be done via the reachProbFunction.
+             *                             Otherwise, the model will be instantiated and checked
+             * 
+             * @return true if an violated point as well as a sat point has been found, i.e., the check result is changed to EXISTSOTH
+             */
+            bool checkPoint(ParameterRegion& region, std::map<VariableType, CoefficientType>const& point, bool viaReachProbFunction=false);
+            
+            /*!
+             * Returns the reachability probability function. 
+             * If it is not yet available, it is computed via state elimination.
+             * After that, the function is available and for the next call of this method it will not be computed again.
+             */
+            ParametricType getReachProbFunction();
+            
+            /*!
+             * Starts the SMTSolver to get the result.
+             * The current regioncheckresult of the region should be EXISTSSAT or EXISTVIOLATED.
+             * Otherwise, a sampingPoint will be computed.
+             * True is returned iff the solver was successful (i.e., it returned sat or unsat)
+             * A Sat- or Violated point is set, if the solver has found one (not yet implemented!).
+             * The region checkResult of the given region is changed accordingly.
+             */
+            bool checkFullSmt(ParameterRegion& region); 
+            
+            //initializes the given solver which can later be used to give an exact result regarding the whole model.
+            void initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, ParametricType const& reachProbFunction, storm::logic::ProbabilityOperatorFormula const& formula);
+            
+            /*!
+             * Returns true iff the given value satisfies the bound given by the specified property
+             */
+            template <typename ValueType>
+            bool valueIsInBoundOfFormula(ValueType value);
+            
+            // The model this model checker is supposed to analyze.
+            storm::models::sparse::Dtmc<ParametricType> const& model;
+            
+            //classes that provide auxilliary functions
+            // Instance of an elimination model checker to access its functions
+            storm::modelchecker::SparseDtmcEliminationModelChecker<ParametricType> eliminationModelChecker;
+            // comparators that can be used to compare constants.
+            storm::utility::ConstantsComparator<ParametricType> parametricTypeComparator;
+            storm::utility::ConstantsComparator<ConstantType> constantTypeComparator;
+            
+            //the following members depend on the currently specified formula:
+            //the currently specified formula
+            std::unique_ptr<storm::logic::ProbabilityOperatorFormula> probabilityOperatorFormula;
+            // the original model after states with constant transitions have been eliminated
+            std::shared_ptr<storm::models::sparse::Dtmc<ParametricType>> simplifiedModel;
+            // the model that can be instantiated to check the value at a certain point
+            std::shared_ptr<SamplingModel> samplingModel;
+            // the model that  is used to approximate the probability values
+            std::shared_ptr<ApproximationModel> approximationModel;
+            // The  function for the reachability probability in the initial state 
+            ParametricType reachProbFunction;
+            // a flag that is true if there are only linear functions at transitions of the model
+            bool hasOnlyLinearFunctions;
+            // a flag that is true iff the reachability probability function has been computed
+            bool isReachProbFunctionComputed;
+            // a flag that is true iff the resulting reachability probability is constant
+            bool isResultConstant;
+            // the smt solver that is used to prove properties with the help of the reachProbFunction
+            std::shared_ptr<storm::solver::Smt2SmtSolver> smtSolver;
+            
+            // runtimes and other information for statistics. 
+            uint_fast64_t numOfCheckedRegions;
+            uint_fast64_t numOfRegionsSolvedThroughSampling;
+            uint_fast64_t numOfRegionsSolvedThroughApproximation;
+            uint_fast64_t numOfRegionsSolvedThroughSubsystemSmt;
+            uint_fast64_t numOfRegionsSolvedThroughFullSmt;
+            uint_fast64_t numOfRegionsExistsBoth;
+            uint_fast64_t numOfRegionsAllSat;
+            uint_fast64_t numOfRegionsAllViolated;
+            
+            std::chrono::high_resolution_clock::duration timePreprocessing;
+            std::chrono::high_resolution_clock::duration timeInitialStateElimination;
+            std::chrono::high_resolution_clock::duration timeComputeReachProbFunction;
+            std::chrono::high_resolution_clock::duration timeCheckRegion;
+            std::chrono::high_resolution_clock::duration timeSampling;
+            std::chrono::high_resolution_clock::duration timeApproximation;
+            std::chrono::high_resolution_clock::duration timeMDPBuild;
+            std::chrono::high_resolution_clock::duration timeSubsystemSmt;
+            std::chrono::high_resolution_clock::duration timeFullSmt;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_ */
diff --git a/src/solver/Smt2SmtSolver.cpp b/src/solver/Smt2SmtSolver.cpp
index 941fe7a3e..2081dd379 100644
--- a/src/solver/Smt2SmtSolver.cpp
+++ b/src/solver/Smt2SmtSolver.cpp
@@ -91,7 +91,7 @@ namespace storm {
         }
         
         void Smt2SmtSolver::add(carl::Constraint<storm::RationalFunction> const& constraint) {
-            add(constraint.lhs(), constraint.rel());
+            add(constraint.lhs(), constraint.relation());
         }
         
         void Smt2SmtSolver::add(carl::Constraint<storm::RawPolynomial> const& constraint) {
@@ -146,7 +146,7 @@ namespace storm {
             if (processIdOfSolver!=0){
                 auto solverOutput = readSolverOutput();
                 STORM_LOG_THROW(solverOutput.size()==1, storm::exceptions::UnexpectedException, "expected a single line of output after smt2 command ( check-sat ). Got " + std::to_string(solverOutput.size()) + " lines of output instead.");
-                solverOutput[0].erase(std::remove_if(solverOutput[0].begin(), solverOutput[0].end(), isspace), solverOutput[0].end()); //remove spaces
+                solverOutput[0].erase(std::remove_if(solverOutput[0].begin(), solverOutput[0].end(), ::isspace), solverOutput[0].end()); //remove spaces
                 if(solverOutput[0]=="sat") return SmtSolver::CheckResult::Sat;
                 if(solverOutput[0]=="unsat") return SmtSolver::CheckResult::Unsat;
                 if(solverOutput[0]=="unknown") return SmtSolver::CheckResult::Unknown;
@@ -262,7 +262,7 @@ namespace storm {
                 if(expectSuccess){
                     auto output=readSolverOutput();
                     STORM_LOG_THROW(output.size()==1, storm::exceptions::UnexpectedException, "expected a single success response after smt2 command " + smt2Command + ". Got " + std::to_string(output.size()) + " lines of output instead.");
-                    output[0].erase(std::remove_if(output[0].begin(), output[0].end(), isspace), output[0].end());
+                    output[0].erase(std::remove_if(output[0].begin(), output[0].end(), ::isspace), output[0].end());
                     STORM_LOG_THROW(output[0]=="success", storm::exceptions::UnexpectedException, "expected <<success>> response after smt2 command " + smt2Command + ". Got <<" + output[0] + ">> instead");
                 }
             }
diff --git a/src/utility/cli.h b/src/utility/cli.h
index ebeb28dd2..41b1e6670 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -66,8 +66,8 @@ log4cplus::Logger printer;
 // Headers for model checking.
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
-#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
-#include "src/utility/regions.h"
+#include "src/modelchecker/region/ParameterRegion.h"
+#include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
@@ -490,7 +490,7 @@ namespace storm {
                 if(settings.isParametricRegionSet()){
                     std::cout << std::endl << "Model checking property: " << *formula << " for all parameters in the given regions." << std::endl;
                     
-                    auto regions=storm::utility::regions::RegionParser<storm::RationalFunction, double>::getRegionsFromSettings();                    
+                    auto regions=storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::ParameterRegion::getRegionsFromSettings();                    
                     storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction, double> modelchecker(*dtmc);
                     if (modelchecker.canHandle(*formula.get())) {
                         modelchecker.specifyFormula(*formula.get());
@@ -499,10 +499,9 @@ namespace storm {
                     else {
                         STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The parametric region check engine currently does not support this property.");
                     }
-                    std::cout << "... done." << std::endl;
-                    for(auto const& reg : regions){
-                        std::cout << reg.toString() << "      Result: " << reg.checkResultToString() << std::endl;
-                    }
+                    //for(auto const& reg : regions){
+                    //    std::cout << reg.toString() << "      Result: " << reg.checkResultToString() << std::endl;
+                    //}
                     modelchecker.printStatisticsToStream(std::cout);
                     
                 }
diff --git a/src/utility/regions.cpp b/src/utility/regions.cpp
index 7fffea186..47af0dfbb 100644
--- a/src/utility/regions.cpp
+++ b/src/utility/regions.cpp
@@ -8,99 +8,26 @@
 #include <string>
 
 #include "src/utility/regions.h"
+#include "src/utility/constants.h"
+#include "src/utility/macros.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/solver/Smt2SmtSolver.h"
+#include "src/exceptions/IllegalArgumentException.h"
 #include "src/exceptions/NotImplementedException.h"
-#include "src/exceptions/InvalidArgumentException.h"
-#include "adapters/CarlAdapter.h"
-#include "exceptions/InvalidSettingsException.h"
-#include "parser/MappedFile.h"
+
+
 
 namespace storm {
     namespace utility{
         namespace regions {
             
-            template<typename ParametricType, typename ConstantType>
-            void RegionParser<ParametricType, ConstantType>::parseParameterBounds(
-                    std::map<VariableType, CoefficientType>& lowerBounds,
-                    std::map<VariableType, CoefficientType>& upperBounds,
-                    std::string const& parameterBoundsString,
-                    double const precision){
-                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
-                
-                std::string::size_type positionOfFirstRelation = parameterBoundsString.find("<=");
-                STORM_LOG_THROW(positionOfFirstRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find  a '<=' after the first number");
-                std::string::size_type positionOfSecondRelation = parameterBoundsString.find("<=", positionOfFirstRelation+2);
-                STORM_LOG_THROW(positionOfSecondRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find  a '<=' after the parameter");
-                
-                std::string parameter=parameterBoundsString.substr(positionOfFirstRelation+2,positionOfSecondRelation-(positionOfFirstRelation+2));
-                //removes all whitespaces from the parameter string:
-                parameter.erase(std::remove_if(parameter.begin(), parameter.end(), isspace), parameter.end());
-                STORM_LOG_THROW(parameter.length()>0, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundsString << " I could not find a parameter");
-                double lowerBound, upperBound;
-                try{
-                    lowerBound=std::stod(parameterBoundsString.substr(0,positionOfFirstRelation));
-                    upperBound=std::stod(parameterBoundsString.substr(positionOfSecondRelation+2));
-                }
-                catch (std::exception const& exception) {
-                    STORM_LOG_ERROR("Failed to parse the region: " << parameterBoundsString << ". The correct format for regions is lowerBound<=parameter<=upperbound");
-                    throw exception;
-                }
-                
-                VariableType var = getVariableFromString<VariableType>(parameter);
-                CoefficientType lb = convertNumber<double, CoefficientType>(lowerBound, true, actualPrecision);
-                STORM_LOG_WARN_COND((lowerBound==convertNumber<CoefficientType, double>(lb, true, actualPrecision)), "The lower bound of '"<< parameterBoundsString << "' could not be parsed accurately. Increase precision?");
-                CoefficientType ub = convertNumber<double, CoefficientType>(upperBound, false, actualPrecision);
-                STORM_LOG_WARN_COND((upperBound==convertNumber<CoefficientType, double>(ub, true, actualPrecision)), "The upper bound of '"<< parameterBoundsString << "' could not be parsed accurately. Increase precision?");
-                lowerBounds.emplace(std::make_pair(var, lb));  
-                upperBounds.emplace(std::make_pair(var, ub));
-               // std::cout << "parsed bounds " << parameterBoundsString << ": lb=" << lowerBound << " ub=" << upperBound << " param='" << parameter << "' precision=" << actualPrecision << std::endl;
-            }
-            
-            template<typename ParametricType, typename ConstantType>
-            typename RegionParser<ParametricType, ConstantType>::ParameterRegion RegionParser<ParametricType, ConstantType>::parseRegion(std::string const& regionString, double precision){
-                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
-                std::map<VariableType, CoefficientType> lowerBounds;
-                std::map<VariableType, CoefficientType> upperBounds;
-                std::vector<std::string> parameterBounds;
-                boost::split(parameterBounds, regionString, boost::is_any_of(","));
-                for(auto const& parameterBound : parameterBounds){
-                    RegionParser<ParametricType, ConstantType>::parseParameterBounds(lowerBounds, upperBounds, parameterBound, actualPrecision);
-                }
-                return ParameterRegion(lowerBounds, upperBounds);
-            }
-            
-            template<typename ParametricType, typename ConstantType>
-            std::vector<typename RegionParser<ParametricType, ConstantType>::ParameterRegion> RegionParser<ParametricType, ConstantType>::parseMultipleRegions(std::string const& regionsString, double precision){
-                double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
-                std::vector<ParameterRegion> result;
-                std::vector<std::string> regionsStrVec;
-                boost::split(regionsStrVec, regionsString, boost::is_any_of(";"));
-                for(auto const& regionStr : regionsStrVec){
-                    if(!std::all_of(regionStr.begin(),regionStr.end(),isspace)){ //skip this string if it only consists of space
-                    result.emplace_back(RegionParser<ParametricType, ConstantType>::parseRegion(regionStr, actualPrecision));
-                    }
-                }
-                return result;
-            }
-            
-            template<typename ParametricType, typename ConstantType>
-            std::vector<typename RegionParser<ParametricType, ConstantType>::ParameterRegion> RegionParser<ParametricType, ConstantType>::getRegionsFromSettings(double precision){
-                STORM_LOG_THROW(storm::settings::regionSettings().isRegionsSet() || storm::settings::regionSettings().isRegionFileSet(), storm::exceptions::InvalidSettingsException, "Tried to obtain regions from the settings but no regions are specified.");
-                STORM_LOG_THROW(!(storm::settings::regionSettings().isRegionsSet() && storm::settings::regionSettings().isRegionFileSet()), storm::exceptions::InvalidSettingsException, "Regions are specified via file AND cmd line. Only one option is allowed.");
-                
-                std::string regionsString;
-                if(storm::settings::regionSettings().isRegionsSet()){
-                    regionsString = storm::settings::regionSettings().getRegionsFromCmdLine();
-                }
-                else{
-                    //if we reach this point we can assume that the region is given as a file.
-                    STORM_LOG_THROW(storm::parser::MappedFile::fileExistsAndIsReadable(storm::settings::regionSettings().getRegionFilePath().c_str()), storm::exceptions::InvalidSettingsException, "The path to the file in which the regions are specified is not valid.");
-                    storm::parser::MappedFile mf(storm::settings::regionSettings().getRegionFilePath().c_str());
-                    regionsString = std::string(mf.getData(),mf.getDataSize());
-                }
-                return RegionParser<ParametricType, ConstantType>::parseMultipleRegions(regionsString,precision);
+            template<>
+            double convertNumber<double, double>(double const& number, bool const& roundDown, double const& precision){
+                return number;
             }
             
-            
+#ifdef STORM_HAVE_CARL
             template<>
             storm::Coefficient convertNumber<double, storm::Coefficient>(double const& number, bool const& roundDown, double const& precision){
                 double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision);
@@ -126,10 +53,6 @@ namespace storm {
                 return cln::double_approx(number);
             }
             
-            template<>
-            double convertNumber<double, double>(double const& number, bool const& roundDown, double const& precision){
-                return number;
-            }
             
             template<>
             cln::cl_RA convertNumber<cln::cl_RA, cln::cl_RA>(cln::cl_RA const& number, bool const& roundDown, double const& precision){
@@ -140,7 +63,7 @@ namespace storm {
             template<>
             storm::Variable getVariableFromString<storm::Variable>(std::string variableString){
                 storm::Variable const& var = carl::VariablePool::getInstance().findVariableWithName(variableString);
-                STORM_LOG_THROW(var!=carl::Variable::NO_VARIABLE, storm::exceptions::InvalidArgumentException, "Variable '" + variableString + "' could not be found.");
+                STORM_LOG_THROW(var!=carl::Variable::NO_VARIABLE, storm::exceptions::IllegalArgumentException, "Variable '" + variableString + "' could not be found.");
                 return var;
             }
             
@@ -162,9 +85,8 @@ namespace storm {
                 }
                 
                 storm::Variable const& var = carl::VariablePool::getInstance().findVariableWithName(variableName);
-                //STORM_LOG_THROW(var==carl::Variable::NO_VARIABLE, storm::exceptions::InvalidArgumentException, "Tried to create a new variable but the name " << variableName << " is already in use.");
                 if(var!=carl::Variable::NO_VARIABLE){
-                    STORM_LOG_THROW(var.getType()==carlVarType, storm::exceptions::InvalidArgumentException, "Tried to create a new variable but the name " << variableName << " is already in use for a variable of a different sort.");
+                    STORM_LOG_THROW(var.getType()==carlVarType, storm::exceptions::IllegalArgumentException, "Tried to create a new variable but the name " << variableName << " is already in use for a variable of a different sort.");
                     return var;
                 }
                 
@@ -177,15 +99,12 @@ namespace storm {
             }
                         
             template<>
-            typename storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::CoefficientType evaluateFunction<storm::RationalFunction, double>(
-                    storm::RationalFunction const& function, 
-                    std::map<typename storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::VariableType,
-                             typename storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::CoefficientType> const& point){
+            CoefficientType<storm::RationalFunction> evaluateFunction<storm::RationalFunction>(storm::RationalFunction const& function, std::map<VariableType<storm::RationalFunction>, CoefficientType<storm::RationalFunction>> const& point){
                 return function.evaluate(point);
             }
             
             template<>
-            typename storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::CoefficientType getConstantPart<storm::RationalFunction, double>(storm::RationalFunction const& function){
+            CoefficientType<storm::RationalFunction> getConstantPart<storm::RationalFunction>(storm::RationalFunction const& function){
                 return function.constantPart();
             }
             
@@ -201,7 +120,7 @@ namespace storm {
             }
             
             template<>
-            void gatherOccurringVariables<storm::RationalFunction, storm::Variable>(storm::RationalFunction const& function, std::set<storm::Variable>& variableSet){
+            void gatherOccurringVariables<storm::RationalFunction>(storm::RationalFunction const& function, std::set<VariableType<storm::RationalFunction>>& variableSet){
                 function.gatherVariables(variableSet);
             }
             
@@ -211,19 +130,19 @@ namespace storm {
                 storm::CompareRelation compRel;
                 switch (relation){
                     case storm::logic::ComparisonType::Greater:
-                        compRel=storm::CompareRelation::GT;
+                        compRel=storm::CompareRelation::GREATER;
                         break;
                     case storm::logic::ComparisonType::GreaterEqual:
                         compRel=storm::CompareRelation::GEQ;
                         break;
                     case storm::logic::ComparisonType::Less:
-                        compRel=storm::CompareRelation::LT;
+                        compRel=storm::CompareRelation::LESS;
                         break;
                     case storm::logic::ComparisonType::LessEqual:
                         compRel=storm::CompareRelation::LEQ;
                         break;
                     default:
-                        STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
+                        STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "the comparison relation of the formula is not supported");
                 }        
                //Note: this only works if numerators and denominators are positive...
                 carl::Constraint<storm::Polynomial> constraint((leftHandSide.nominator() * rightHandSide.denominator()) - (rightHandSide.nominator() * leftHandSide.denominator()), compRel);
@@ -235,19 +154,19 @@ namespace storm {
                 storm::CompareRelation compRel;
                 switch (relation){
                     case storm::logic::ComparisonType::Greater:
-                        compRel=storm::CompareRelation::GT;
+                        compRel=storm::CompareRelation::GREATER;
                         break;
                     case storm::logic::ComparisonType::GreaterEqual:
                         compRel=storm::CompareRelation::GEQ;
                         break;
                     case storm::logic::ComparisonType::Less:
-                        compRel=storm::CompareRelation::LT;
+                        compRel=storm::CompareRelation::LESS;
                         break;
                     case storm::logic::ComparisonType::LessEqual:
                         compRel=storm::CompareRelation::LEQ;
                         break;
                     default:
-                        STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported");
+                        STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "the comparison relation of the formula is not supported");
                 }
                 storm::RawPolynomial leftHandSide(variable);
                 leftHandSide -= bound;
@@ -271,33 +190,7 @@ namespace storm {
                 std::shared_ptr<carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>> cache(new carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>());
                 return storm::RationalFunction(storm::RationalFunction::PolyType(storm::RationalFunction::PolyType::PolyType(initialValue), cache));
             }
-
-            
-            //explicit instantiations
-       template double convertNumber<double, double>(double const& number, bool const& roundDown, double const& precision);
-       
-#ifdef STORM_HAVE_CARL
-       template class RegionParser<storm::RationalFunction, double>;
-       
-       template storm::RationalFunction convertNumber<double, storm::RationalFunction>(double const& number, bool const& roundDown, double const& precision);
-       template storm::Coefficient convertNumber<double, storm::Coefficient>(double const& number, bool const& roundDown, double const& precision);
-       template double convertNumber<cln::cl_RA, double>(storm::Coefficient const& number, bool const& roundDown, double const& precision);
-       template cln::cl_RA convertNumber<cln::cl_RA, cln::cl_RA>(cln::cl_RA const& number, bool const& roundDown, double const& precision);
-       
-       template storm::Variable getVariableFromString<storm::Variable>(std::string variableString);
-       template std::string getVariableName<storm::Variable>(storm::Variable variable);
-       
-       template bool functionIsLinear<storm::RationalFunction>(storm::RationalFunction const& function);
-       
-       template void gatherOccurringVariables<storm::RationalFunction, storm::Variable>(storm::RationalFunction const& function, std::set<storm::Variable>& variableSet);
-       
-       template void addGuardedConstraintToSmtSolver<storm::solver::Smt2SmtSolver, storm::RationalFunction, storm::Variable>(std::shared_ptr<storm::solver::Smt2SmtSolver> solver,storm::Variable const& guard, storm::RationalFunction const& leftHandSide, storm::logic::ComparisonType relation, storm::RationalFunction const& rightHandSide);
-       template void addParameterBoundsToSmtSolver<storm::solver::Smt2SmtSolver, storm::Variable, cln::cl_RA>(std::shared_ptr<storm::solver::Smt2SmtSolver> solver, storm::Variable const& variable, storm::logic::ComparisonType relation, cln::cl_RA const& bound);
-       template void addBoolVariableToSmtSolver<storm::solver::Smt2SmtSolver, storm::Variable>(std::shared_ptr<storm::solver::Smt2SmtSolver> solver, storm::Variable const& variable, bool value);
-       
-       template storm::RationalFunction getNewFunction<storm::RationalFunction, storm::Coefficient>(storm::Coefficient initialValue);
-       template storm::RationalFunction getNewFunction<storm::RationalFunction, storm::Variable>(storm::Variable initialValue);
-#endif 
+#endif
         }
     }
 }
diff --git a/src/utility/regions.h b/src/utility/regions.h
index e9e338c3d..b20eb098e 100644
--- a/src/utility/regions.h
+++ b/src/utility/regions.h
@@ -12,8 +12,13 @@
 #ifndef STORM_UTILITY_REGIONS_H
 #define	STORM_UTILITY_REGIONS_H
 
-#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
+#include <type_traits>
+#include <map>
+#include <set>
+#include <memory>
+
 #include "src/logic/ComparisonType.h"
+#include <src/adapters/CarlAdapter.h>
 
 // Forward-declare region modelchecker  class.
     namespace storm {
@@ -26,62 +31,22 @@
 namespace storm {
     namespace utility{
         namespace regions {
-            template<typename ParametricType, typename ConstantType>
-            class RegionParser{
-                public:
-                    typedef typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::ParameterRegion ParameterRegion;
-                    typedef typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::VariableType VariableType;
-                    typedef typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType CoefficientType;
-                    
-                /*
-                 * Can be used to parse a single parameter with its bounds from a string of the form "0.3<=p<=0.5".
-                 * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
-                 * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
-                 * If no precision is given, the one from the settings is used.
-                 * The results will be inserted in the given maps
-                 * 
-                 */
-                static void parseParameterBounds( 
-                        std::map<VariableType, CoefficientType>& lowerBounds,
-                        std::map<VariableType, CoefficientType>& upperBounds,
-                        std::string const& parameterBoundsString,
-                        double const precision=0.0
-                );
-
-                /*
-                 * Can be used to parse a single region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7".
-                 * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
-                 * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
-                 * If no precision is given, the one from the settings is used.
-                 * 
-                 */
-                static ParameterRegion parseRegion(
-                        std::string const& regionString,
-                        double precision=0.0);
-
-                /*
-                 * Can be used to parse a vector of region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7;0.1<=p<=0.3,0.2<=q<=0.4".
-                 * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
-                 * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
-                 * If no precision is given, the one from the settings is used.
-                 * 
-                 */
-                static std::vector<ParameterRegion> parseMultipleRegions(
-                        std::string const& regionsString,
-                        double precision=0.0);
-            
-
-                /*
-                 * Retrieves the regions that are specified in the settings.
-                 * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType.
-                 * According to the given precision, the lower bound may be rounded down and the upper bound may be rounded up.
-                 * If no precision is given, the one from the settings is used.
-                 * 
-                 */
-                static std::vector<ParameterRegion> getRegionsFromSettings(double precision=0.0);
             
-            };
+#ifdef STORM_HAVE_CARL
+            template<typename FunctionType>
+            using VariableType    = typename std::conditional<(std::is_same<FunctionType, storm::RationalFunction>::value), storm::Variable, std::nullptr_t>::type;
+            template<typename FunctionType>
+            using CoefficientType = typename std::conditional<(std::is_same<FunctionType, storm::RationalFunction>::value), storm::Coefficient, std::nullptr_t>::type;
+#else
+            template<typename Functiontype>
+            using VariableType = std::nullptr_t;
+            template<typename Functiontype>
+            using CoefficientType = std::nullptr_t;
+#endif
+       //     template<typename FunctionType>
+        //    using PointType = std::map<VariableType<FunctionType>, CoefficientType<FunctionType>>;
             
+            enum class VariableSort {VS_BOOL, VS_REAL, VS_INT};
             
             /*
              * Converts a number from one type to a number from the other.
@@ -94,50 +59,47 @@ namespace storm {
              * retrieves the variable object from the given string
              * Throws an exception if variable not found
              */
-            template<typename VariableType>
-            VariableType getVariableFromString(std::string variableString);
-            
-            
-            enum class VariableSort {VS_BOOL, VS_REAL, VS_INT};
+            template<typename VarType>
+            VarType getVariableFromString(std::string variableString);
             
             /*
              * Creates a new variable with the given name and the given sort
              * If there is already a variable with that name, that variable is returned.
              * An exception is thrown if there already is a variable with the given name, but with a different sort.
              */
-            template<typename VariableType>
-            VariableType getNewVariable(std::string variableName, VariableSort sort);
+            template<typename VarType>
+            VarType getNewVariable(std::string variableName, VariableSort sort);
             
             /*
              * retrieves the variable name  from the given variable
              */
-            template<typename VariableType>
-            std::string getVariableName(VariableType variable);
+            template<typename VarType>
+            std::string getVariableName(VarType variable);
             
             /*
              * evaluates the given function at the given point and returns the result
              */
-            template<typename ParametricType, typename ConstantType>
-            typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType evaluateFunction(ParametricType const& function, std::map<typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::VariableType, typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType> const& point);
+            template<typename FunctionType>
+            CoefficientType<FunctionType> evaluateFunction(FunctionType const& function, std::map<VariableType<FunctionType>, CoefficientType<FunctionType>> const& point);
             
             /*!
              * retrieves the constant part of the given function.
              * If the function is constant, then the result is the same value as the original function
              */
-            template<typename ParametricType, typename ConstantType>
-            typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType getConstantPart(ParametricType const& function);
+            template<typename FunctionType>
+            CoefficientType<FunctionType> getConstantPart(FunctionType const& function);
             
             /*!
              * Returns true if the function is rational. Note that the function might be simplified.
              */
-            template<typename ParametricType>
-            bool functionIsLinear(ParametricType const& function);
+            template<typename FunctionType>
+            bool functionIsLinear(FunctionType const& function);
             
             /*!
              *  Add all variables that occur in the given function to the the given set
              */
-            template<typename ParametricType, typename VariableType>
-            void gatherOccurringVariables(ParametricType const& function, std::set<VariableType>& variableSet);
+            template<typename FunctionType>
+            void gatherOccurringVariables(FunctionType const& function, std::set<VariableType<FunctionType>>& variableSet);
             
             
             /*!
@@ -150,29 +112,28 @@ namespace storm {
              * @param relation relation of the constraint
              * @param rightHandSide right hand side of the constraint
              */
-            template<typename SolverType, typename ParametricType, typename VariableType>
-            void addGuardedConstraintToSmtSolver(std::shared_ptr<SolverType> solver, VariableType const& guard, ParametricType const& leftHandSide, storm::logic::ComparisonType relation, ParametricType const& rightHandSide);
+            template<typename SolverType, typename FunctionType, typename VarType>
+            void addGuardedConstraintToSmtSolver(std::shared_ptr<SolverType> solver, VarType const& guard, FunctionType const& leftHandSide, storm::logic::ComparisonType relation, FunctionType const& rightHandSide);
             
             /*!
              * Adds the given constraint to the given Smt solver
              * The constraint is of the form 'variable relation bound'
              */
-            template<typename SolverType, typename VariableType, typename BoundType>
-            void addParameterBoundsToSmtSolver(std::shared_ptr<SolverType> solver, VariableType const& variable, storm::logic::ComparisonType relation, BoundType const& bound);
+            template<typename SolverType, typename VarType, typename BoundType>
+            void addParameterBoundsToSmtSolver(std::shared_ptr<SolverType> solver, VarType const& variable, storm::logic::ComparisonType relation, BoundType const& bound);
             
             /*!
              * Adds the given (boolean) variable to the solver. Can be used to assert that guards are true/false
              */
-            template<typename SolverType, typename VariableType>
-            void addBoolVariableToSmtSolver(std::shared_ptr<SolverType> solver, VariableType const& variable, bool value);
+            template<typename SolverType, typename VarType>
+            void addBoolVariableToSmtSolver(std::shared_ptr<SolverType> solver, VarType const& variable, bool value);
             
             
             /*!
              * Returns a new function that initially has the given value (which might be a constant or a variable)
-             * 
              */
-            template<typename ParametricType, typename ValueType>
-            ParametricType getNewFunction(ValueType initialValue);
+            template<typename FunctionType, typename ValueType>
+            FunctionType getNewFunction(ValueType initialValue);
         }
     }
 }