diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
index 428584385..413f34ac8 100644
--- a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
@@ -1,19 +1,9 @@
 #include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
-#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h"
-#include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
-
-#include "src/storage/dd/CuddOdd.h"
 
-#include "src/utility/macros.h"
-#include "src/utility/graph.h"
+#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h"
+#include "src/modelchecker/csl/helper/HybridCtmcCslHelper.h"
 
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
-#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
-#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
-#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
-
-#include "src/exceptions/InvalidStateException.h"
-#include "src/exceptions/InvalidPropertyException.h"
 
 namespace storm {
     namespace modelchecker {
@@ -31,12 +21,7 @@ namespace storm {
         bool HybridCtmcCslModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const {
             return formula.isCslStateFormula() || formula.isCslPathFormula() || formula.isRewardPathFormula();
         }
-        
-        template<storm::dd::DdType DdType, class ValueType>
-        storm::dd::Add<DdType> HybridCtmcCslModelChecker<DdType, ValueType>::computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector) {
-            return rateMatrix / exitRateVector;
-        }
-        
+                
         template<storm::dd::DdType DdType, class ValueType>
         std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
@@ -44,14 +29,15 @@ namespace storm {
             SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             
-            return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
         }
         
         template<storm::dd::DdType DdType, class ValueType>
         std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), HybridDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector())));
+            
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeNextProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), subResult.getTruthValuesVector());
         }
 
         template<storm::dd::DdType DdType, class ValueType>
@@ -60,34 +46,11 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, class ValueType>
-        storm::dd::Add<DdType> HybridCtmcCslModelChecker<DdType, ValueType>::computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate) {
-            STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << ".");
-            STORM_LOG_DEBUG("Keeping " << maybeStates.getNonZeroCount() << " rows.");
-            
-            // Cut all non-maybe rows/columns from the transition matrix.
-            storm::dd::Add<DdType> uniformizedMatrix = transitionMatrix * maybeStates.toAdd() * maybeStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
-            
-            // Now perform the uniformization.
-            uniformizedMatrix = uniformizedMatrix / model.getManager().getConstant(uniformizationRate);
-            storm::dd::Add<DdType> diagonal = model.getRowColumnIdentity() * maybeStates.toAdd();
-            storm::dd::Add<DdType> diagonalOffset = diagonal;
-            diagonalOffset -= diagonal * (exitRateVector / model.getManager().getConstant(uniformizationRate));
-            uniformizedMatrix += diagonalOffset;
-            
-            return uniformizedMatrix;
-        }
-        
-        template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
 
-            boost::optional<storm::dd::Add<DdType>> modifiedStateRewardVector;
-            if (this->getModel().hasStateRewards()) {
-                modifiedStateRewardVector = this->getModel().getStateRewardVector() / this->getModel().getExitRateVector();
-            }
-            
-            return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative);
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeReachabilityRewards(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
         }
         
         template<storm::dd::DdType DdType, class ValueType>
@@ -106,255 +69,17 @@ namespace storm {
                 upperBound = pathFormula.getDiscreteTimeBound();
             }
             
-            return this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound);
-        }
-        
-        template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, storm::dd::Add<DdType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const {
-            // If the time bounds are [0, inf], we rather call untimed reachability.
-            storm::utility::ConstantsComparator<ValueType> comparator;
-            if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
-                return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
-            }
-            
-            // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0
-            // or t' != inf.
-            
-            // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the
-            // further computations.
-            storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), this->getModel().getTransitionMatrix().notZero(), phiStates, psiStates);
-            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0.getNonZeroCount() << " states with probability greater 0.");
-            storm::dd::Bdd<DdType> statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 && !psiStates;
-            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNonZeroCount() << " 'maybe' states.");
-            
-            if (!statesWithProbabilityGreater0NonPsi.isZero()) {
-                if (comparator.isZero(upperBound)) {
-                    // In this case, the interval is of the form [0, 0].
-                    return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), psiStates.toAdd()));
-                } else {
-                    if (comparator.isZero(lowerBound)) {
-                        // In this case, the interval is of the form [0, t].
-                        // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier.
-                        
-                        // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
-                        ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRates).getMax();
-                        STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                        
-                        // Compute the uniformized matrix.
-                        storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates, statesWithProbabilityGreater0NonPsi, uniformizationRate);
-                        
-                        // Compute the vector that is to be added as a compensation for removing the absorbing states.
-                        storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * this->getModel().getTransitionMatrix() * psiStates.swapVariables(this->getModel().getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(this->getModel().getColumnVariables()) / this->getModel().getManager().getConstant(uniformizationRate);
-                        
-                        // Create an ODD for the translation to an explicit representation.
-                        storm::dd::Odd<DdType> odd(statesWithProbabilityGreater0NonPsi);
-                        
-                        // Convert the symbolic parts to their explicit representation.
-                        storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                        std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
-                        
-                        // Finally compute the transient probabilities.
-                        std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
-                        std::vector<ValueType> subresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, *this->linearEquationSolverFactory);
-
-                        return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(),
-                                                                                                      (psiStates || !statesWithProbabilityGreater0) && this->getModel().getReachableStates(),
-                                                                                                      psiStates.toAdd(),
-                                                                                                      statesWithProbabilityGreater0NonPsi,
-                                                                                                      odd, subresult));
-                    } else if (comparator.isInfinity(upperBound)) {
-                        // In this case, the interval is of the form [t, inf] with t != 0.
-                        
-                        // Start by computing the (unbounded) reachability probabilities of reaching psi states while
-                        // staying in phi states.
-                        std::unique_ptr<CheckResult> unboundedResult = HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
-                        
-                        // Compute the set of relevant states.
-                        storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
-
-                        // Filter the unbounded result such that it only contains values for the relevant states.
-                        unboundedResult->filter(SymbolicQualitativeCheckResult<DdType>(this->getModel().getReachableStates(), relevantStates));
-
-                        // Build an ODD for the relevant states.
-                        storm::dd::Odd<DdType> odd(relevantStates);
-                        
-                        std::vector<ValueType> result;
-                        if (unboundedResult->isHybridQuantitativeCheckResult()) {
-                            std::unique_ptr<CheckResult> explicitUnboundedResult = unboundedResult->asHybridQuantitativeCheckResult<DdType>().toExplicitQuantitativeCheckResult();
-                            result = std::move(explicitUnboundedResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
-                        } else {
-                            STORM_LOG_THROW(unboundedResult->isSymbolicQuantitativeCheckResult(), storm::exceptions::InvalidStateException, "Expected check result of different type.");
-                            result = unboundedResult->asSymbolicQuantitativeCheckResult<DdType>().getValueVector().template toVector<ValueType>(odd);
-                        }
-                        
-                        // Determine the uniformization rate for the transient probability computation.
-                        ValueType uniformizationRate = 1.02 * (relevantStates.toAdd() * exitRates).getMax();
-                        
-                        // Compute the uniformized matrix.
-                        storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  relevantStates, uniformizationRate);
-                        storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                        
-                        // Compute the transient probabilities.
-                        result = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, *this->linearEquationSolverFactory);
-                        
-                        return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, result));
-                    } else {
-                        // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
-                        
-                        if (lowerBound != upperBound) {
-                            // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'.
-
-                            // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
-                            ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRates).getMax();
-                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                        
-                            // Compute the (first) uniformized matrix.
-                            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0NonPsi, uniformizationRate);
-                            
-                            // Create the one-step vector.
-                            storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * this->getModel().getTransitionMatrix() * psiStates.swapVariables(this->getModel().getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(this->getModel().getColumnVariables()) / this->getModel().getManager().getConstant(uniformizationRate);
-                            
-                            // Build an ODD for the relevant states and translate the symbolic parts to their explicit representation.
-                            storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0NonPsi);
-                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                            std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
-                            
-                            // Compute the transient probabilities.
-                            std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
-                            std::vector<ValueType> subResult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
-
-                            // Transform the explicit result to a hybrid check result, so we can easily convert it to
-                            // a symbolic qualitative format.
-                            HybridQuantitativeCheckResult<DdType> hybridResult(this->getModel().getReachableStates(), psiStates || (!statesWithProbabilityGreater0 && this->getModel().getReachableStates()),
-                                                                               psiStates.toAdd(), statesWithProbabilityGreater0NonPsi, odd, subResult);
-                            
-                            // Compute the set of relevant states.
-                            storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
-                            
-                            // Filter the unbounded result such that it only contains values for the relevant states.
-                            hybridResult.filter(SymbolicQualitativeCheckResult<DdType>(this->getModel().getReachableStates(), relevantStates));
-                            
-                            // Build an ODD for the relevant states.
-                            odd = storm::dd::Odd<DdType>(relevantStates);
-                            
-                            std::unique_ptr<CheckResult> explicitResult = hybridResult.toExplicitQuantitativeCheckResult();
-                            std::vector<ValueType> newSubresult = std::move(explicitResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
-                            
-                            // Then compute the transient probabilities of being in such a state after t time units. For this,
-                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
-                            uniformizationRate =  1.02 * (relevantStates.toAdd() * exitRates).getMax();
-                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                            
-                            // If the lower and upper bounds coincide, we have only determined the relevant states at this
-                            // point, but we still need to construct the starting vector.
-                            if (lowerBound == upperBound) {
-                                odd = storm::dd::Odd<DdType>(relevantStates);
-                                newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
-                            }
-                            
-                            // Finally, we compute the second set of transient probabilities.
-                            uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  relevantStates, uniformizationRate);
-                            explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                            
-                            newSubresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
-                            
-                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, newSubresult));
-                        } else {
-                            // In this case, the interval is of the form [t, t] with t != 0, t != inf.
-
-                            // Build an ODD for the relevant states.
-                            storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0);
-                            
-                            std::vector<ValueType> newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
-                            
-                            // Then compute the transient probabilities of being in such a state after t time units. For this,
-                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
-                            ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0.toAdd() * exitRates).getMax();
-                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                            
-                            // Finally, we compute the second set of transient probabilities.
-                            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0, uniformizationRate);
-                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                            
-                            newSubresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
-                            
-                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !statesWithProbabilityGreater0 && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), statesWithProbabilityGreater0, odd, newSubresult));
-                        }
-                    }
-                }
-            } else {
-                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), psiStates.toAdd()));
-            }
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeBoundedUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, lowerBound, upperBound, *linearEquationSolverFactory);
         }
         
         template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            return this->computeInstantaneousRewardsHelper(rewardPathFormula.getContinuousTimeBound());
-        }
-        
-        template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(double timeBound) const {
-            // Only compute the result if the model has a state-based reward this->getModel().
-            STORM_LOG_THROW(this->getModel().hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Create ODD for the translation.
-            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
-            
-            // Initialize result to state rewards of the this->getModel().
-            std::vector<ValueType> result = this->getModel().getStateRewardVector().template toVector<ValueType>(odd);
-            
-            // If the time-bound is not zero, we need to perform a transient analysis.
-            if (timeBound > 0) {
-                ValueType uniformizationRate = 1.02 * this->getModel().getExitRateVector().getMax();
-                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                
-                storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(),  this->getModel().getReachableStates(), uniformizationRate);
-                
-                storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                result = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory);
-            }
-            
-            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), odd, result));
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeInstantaneousRewards(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardPathFormula.getContinuousTimeBound(), *linearEquationSolverFactory);
         }
 
         template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            return this->computeCumulativeRewardsHelper(rewardPathFormula.getContinuousTimeBound());
-        }
-        
-        template<storm::dd::DdType DdType, class ValueType>
-        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(double timeBound) const {
-            // Only compute the result if the model has a state-based reward this->getModel().
-            STORM_LOG_THROW(this->getModel().hasStateRewards() || this->getModel().hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // If the time bound is zero, the result is the constant zero vector.
-            if (timeBound == 0) {
-                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getAddZero()));
-            }
-            
-            // Otherwise, we need to perform some computations.
-            
-            // Start with the uniformization.
-            ValueType uniformizationRate = 1.02 * this->getModel().getExitRateVector().getMax();
-            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-            
-            // Create ODD for the translation.
-            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
-            
-            // Compute the uniformized matrix.
-            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(),  this->getModel().getReachableStates(), uniformizationRate);
-            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-
-            // Then compute the state reward vector to use in the computation.
-            storm::dd::Add<DdType> totalRewardVector = this->getModel().hasStateRewards() ? this->getModel().getStateRewardVector() : this->getModel().getManager().getAddZero();
-            if (this->getModel().hasTransitionRewards()) {
-                totalRewardVector += (this->getModel().getTransitionMatrix() * this->getModel().getTransitionRewardMatrix()).sumAbstract(this->getModel().getColumnVariables());
-            }
-            std::vector<ValueType> explicitTotalRewardVector = totalRewardVector.template toVector<ValueType>(odd);
-            
-            // Finally, compute the transient probabilities.
-            std::vector<ValueType> result = SparseCtmcCslHelper<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, *this->linearEquationSolverFactory);
-            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result)));
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeCumulativeRewards(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardPathFormula.getContinuousTimeBound(), *linearEquationSolverFactory);
         }
         
         template<storm::dd::DdType DdType, class ValueType>
@@ -362,16 +87,7 @@ namespace storm {
             std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             
-            storm::dd::Add<DdType> probabilityMatrix = this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
-            
-            // Create ODD for the translation.
-            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
-            
-            storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd);
-            std::vector<ValueType> explicitExitRateVector = this->getModel().getExitRateVector().template toVector<ValueType>(odd);
-            
-            std::vector<ValueType> result = SparseCtmcCslHelper<ValueType>::computeLongRunAverage(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), &explicitExitRateVector, qualitative, *this->linearEquationSolverFactory);
-            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result)));
+            return storm::modelchecker::helper::HybridCtmcCslHelper<DdType, ValueType>::computeLongRunAverage(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
         }
         
         // Explicitly instantiate the model checker.
diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.h b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
index b2a94d51e..a1f9cb6d4 100644
--- a/src/modelchecker/csl/HybridCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
@@ -2,9 +2,10 @@
 #define STORM_MODELCHECKER_HYBRIDCTMCCSLMODELCHECKER_H_
 
 #include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+
 #include "src/models/symbolic/Ctmc.h"
+
 #include "src/utility/solver.h"
-#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -20,42 +21,15 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
 
         protected:
             storm::models::symbolic::Ctmc<DdType> const& getModel() const override;
             
         private:
-            /*!
-             * Converts the given rate-matrix into a time-abstract probability matrix.
-             *
-             * @param model The symbolic model.
-             * @param rateMatrix The rate matrix.
-             * @param exitRateVector The exit rate vector of the model.
-             * @return The probability matrix.
-             */
-            static storm::dd::Add<DdType> computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector);
-
-            /*!
-             * Computes the matrix representing the transitions of the uniformized CTMC.
-             *
-             * @param model The symbolic model.
-             * @param transitionMatrix The matrix to uniformize.
-             * @param exitRateVector The exit rate vector.
-             * @param maybeStates The states that need to be considered.
-             * @param uniformizationRate The rate to be used for uniformization.
-             * @return The uniformized matrix.
-             */
-            static storm::dd::Add<DdType> computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate);
-            
-            // The methods that perform the actual checking.
-            std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, storm::dd::Add<DdType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const;
-            std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(double timeBound) const;
-            std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(double timeBound) const;
-            
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
index 99c92ee85..fa7fb98b0 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
@@ -1,208 +1,46 @@
 #include "src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h"
 
-#include <utility>
-#include <vector>
+#include "src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h"
 
-#include "src/storage/StronglyConnectedComponentDecomposition.h"
-
-#include "src/utility/constants.h"
 #include "src/utility/macros.h"
-#include "src/utility/vector.h"
-#include "src/utility/graph.h"
 
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
-#include "src/solver/LpSolver.h"
-
 #include "src/exceptions/InvalidPropertyException.h"
 #include "src/exceptions/NotImplementedException.h"
 
 namespace storm {
     namespace modelchecker {
-        template<typename ValueType>
-        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& MinMaxLinearEquationSolverFactory) : model(model), MinMaxLinearEquationSolverFactory(std::move(MinMaxLinearEquationSolverFactory)) {
+        template<typename SparseMarkovAutomatonModelType>
+        SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::SparseMarkovAutomatonCslModelChecker(SparseMarkovAutomatonModelType const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& minMaxLinearEquationSolverFactory) : SparsePropositionalModelChecker<SparseMarkovAutomatonModelType>(model), minMaxLinearEquationSolverFactory(std::move(minMaxLinearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
-        template<typename ValueType>
-        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model) : model(model), MinMaxLinearEquationSolverFactory(new storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>()) {
+        template<typename SparseMarkovAutomatonModelType>
+        SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::SparseMarkovAutomatonCslModelChecker(SparseMarkovAutomatonModelType const& model) : SparsePropositionalModelChecker<SparseMarkovAutomatonModelType>(model), minMaxLinearEquationSolverFactory(new storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
-        template<typename ValueType>
-        bool SparseMarkovAutomatonCslModelChecker<ValueType>::canHandle(storm::logic::Formula const& formula) const {
+        template<typename SparseMarkovAutomatonModelType>
+        bool SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::canHandle(storm::logic::Formula const& formula) const {
             return formula.isCslStateFormula() || formula.isCslPathFormula() || (formula.isRewardPathFormula() && formula.isReachabilityRewardFormula());
         }
-
-        template<typename ValueType>
-        void SparseMarkovAutomatonCslModelChecker<ValueType>::computeBoundedReachabilityProbabilities(bool min, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory) {
-            // Start by computing four sparse matrices:
-            // * a matrix aMarkovian with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states.
-            // * a matrix aMarkovianToProbabilistic with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states.
-            // * a matrix aProbabilistic with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states.
-            // * a matrix aProbabilisticToMarkovian with all (non-discretized) transitions from probabilistic non-goal states to all Markovian non-goal states.
-            typename storm::storage::SparseMatrix<ValueType> aMarkovian = transitionMatrix.getSubmatrix(true, markovianNonGoalStates, markovianNonGoalStates, true);
-            typename storm::storage::SparseMatrix<ValueType> aMarkovianToProbabilistic = transitionMatrix.getSubmatrix(true, markovianNonGoalStates, probabilisticNonGoalStates);
-            typename storm::storage::SparseMatrix<ValueType> aProbabilistic = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, probabilisticNonGoalStates);
-            typename storm::storage::SparseMatrix<ValueType> aProbabilisticToMarkovian = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, markovianNonGoalStates);
-            
-            // The matrices with transitions from Markovian states need to be digitized.
-            // Digitize aMarkovian. Based on whether the transition is a self-loop or not, we apply the two digitization rules.
-            uint_fast64_t rowIndex = 0;
-            for (auto state : markovianNonGoalStates) {
-                for (auto& element : aMarkovian.getRow(rowIndex)) {
-                    ValueType eTerm = std::exp(-exitRates[state] * delta);
-                    if (element.getColumn() == rowIndex) {
-                        element.setValue((storm::utility::one<ValueType>() - eTerm) * element.getValue() + eTerm);
-                    } else {
-                        element.setValue((storm::utility::one<ValueType>() - eTerm) * element.getValue());
-                    }
-                }
-                ++rowIndex;
-            }
-            
-            // Digitize aMarkovianToProbabilistic. As there are no self-loops in this case, we only need to apply the digitization formula for regular successors.
-            rowIndex = 0;
-            for (auto state : markovianNonGoalStates) {
-                for (auto& element : aMarkovianToProbabilistic.getRow(rowIndex)) {
-                    element.setValue((1 - std::exp(-exitRates[state] * delta)) * element.getValue());
-                }
-                ++rowIndex;
-            }
-            
-            // Initialize the two vectors that hold the variable one-step probabilities to all target states for probabilistic and Markovian (non-goal) states.
-            std::vector<ValueType> bProbabilistic(aProbabilistic.getRowCount());
-            std::vector<ValueType> bMarkovian(markovianNonGoalStates.getNumberOfSetBits());
-            
-            // Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones.
-            std::vector<ValueType> bProbabilisticFixed = transitionMatrix.getConstrainedRowSumVector(probabilisticNonGoalStates, goalStates);
-            std::vector<ValueType> bMarkovianFixed;
-            bMarkovianFixed.reserve(markovianNonGoalStates.getNumberOfSetBits());
-            for (auto state : markovianNonGoalStates) {
-                bMarkovianFixed.push_back(storm::utility::zero<ValueType>());
-                
-                for (auto& element : transitionMatrix.getRowGroup(state)) {
-                    if (goalStates.get(element.getColumn())) {
-                        bMarkovianFixed.back() += (1 - std::exp(-exitRates[state] * delta)) * element.getValue();
-                    }
-                }
-            }
-            
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory.create(aProbabilistic);
-                        
-            // Perform the actual value iteration
-            // * loop until the step bound has been reached
-            // * in the loop:
-            // *    perform value iteration using A_PSwG, v_PS and the vector b where b = (A * 1_G)|PS + A_PStoMS * v_MS
-            //      and 1_G being the characteristic vector for all goal states.
-            // *    perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS
-            std::vector<ValueType> markovianNonGoalValuesSwap(markovianNonGoalValues);
-            std::vector<ValueType> multiplicationResultScratchMemory(aProbabilistic.getRowCount());
-            std::vector<ValueType> aProbabilisticScratchMemory(probabilisticNonGoalValues.size());
-            for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
-                // Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
-                aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
-                storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
-                
-                // Now perform the inner value iteration for probabilistic states.
-                solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
-                
-                // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
-                aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
-                storm::utility::vector::addVectors(bMarkovian, bMarkovianFixed, bMarkovian);
-                aMarkovian.multiplyWithVector(markovianNonGoalValues, markovianNonGoalValuesSwap);
-                std::swap(markovianNonGoalValues, markovianNonGoalValuesSwap);
-                storm::utility::vector::addVectors(markovianNonGoalValues, bMarkovian, markovianNonGoalValues);
-            }
-            
-            // After the loop, perform one more step of the value iteration for PS states.
-            aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
-            storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
-            solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
-        }
         
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair) const {
-            // 'Unpack' the bounds to make them more easily accessible.
-            double lowerBound = boundsPair.first;
-            double upperBound = boundsPair.second;
-            
-            // Get some data fields for convenient access.
-            typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = model.getTransitionMatrix();
-            std::vector<ValueType> const& exitRates = model.getExitRates();
-            storm::storage::BitVector const& markovianStates = model.getMarkovianStates();
-            
-            // (1) Compute the accuracy we need to achieve the required error bound.
-            ValueType maxExitRate = model.getMaximalExitRate();
-            ValueType delta = (2 * storm::settings::generalSettings().getPrecision()) / (upperBound * maxExitRate * maxExitRate);
-            
-            // (2) Compute the number of steps we need to make for the interval.
-            uint_fast64_t numberOfSteps = static_cast<uint_fast64_t>(std::ceil((upperBound - lowerBound) / delta));
-            STORM_LOG_INFO("Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [" << lowerBound << ", " << upperBound << "]." << std::endl);
-            
-            // (3) Compute the non-goal states and initialize two vectors
-            // * vProbabilistic holds the probability values of probabilistic non-goal states.
-            // * vMarkovian holds the probability values of Markovian non-goal states.
-            storm::storage::BitVector const& markovianNonGoalStates = markovianStates & ~psiStates;
-            storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~psiStates;
-            std::vector<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
-            std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
-            
-            computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, psiStates, markovianNonGoalStates, probabilisticNonGoalStates, vMarkovian, vProbabilistic, delta, numberOfSteps, *MinMaxLinearEquationSolverFactory);
-            
-            // (4) If the lower bound of interval was non-zero, we need to take the current values as the starting values for a subsequent value iteration.
-            if (lowerBound != storm::utility::zero<ValueType>()) {
-                std::vector<ValueType> vAllProbabilistic((~markovianStates).getNumberOfSetBits());
-                std::vector<ValueType> vAllMarkovian(markovianStates.getNumberOfSetBits());
-                
-                // Create the starting value vectors for the next value iteration based on the results of the previous one.
-                storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, psiStates % ~markovianStates, storm::utility::one<ValueType>());
-                storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, ~psiStates % ~markovianStates, vProbabilistic);
-                storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, psiStates % markovianStates, storm::utility::one<ValueType>());
-                storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, ~psiStates % markovianStates, vMarkovian);
-                
-                // Compute the number of steps to reach the target interval.
-                numberOfSteps = static_cast<uint_fast64_t>(std::ceil(lowerBound / delta));
-                STORM_LOG_INFO("Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [0, " << lowerBound << "]." << std::endl);
-                
-                // Compute the bounded reachability for interval [0, b-a].
-                computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, storm::storage::BitVector(model.getNumberOfStates()), markovianStates, ~markovianStates, vAllMarkovian, vAllProbabilistic, delta, numberOfSteps, *MinMaxLinearEquationSolverFactory);
-                
-                // Create the result vector out of vAllProbabilistic and vAllMarkovian and return it.
-                std::vector<ValueType> result(model.getNumberOfStates());
-                storm::utility::vector::setVectorValues(result, ~markovianStates, vAllProbabilistic);
-                storm::utility::vector::setVectorValues(result, markovianStates, vAllMarkovian);
-                
-                return result;
-            } else {
-                // Create the result vector out of 1_G, vProbabilistic and vMarkovian and return it.
-                std::vector<ValueType> result(model.getNumberOfStates());
-                storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
-                storm::utility::vector::setVectorValues(result, probabilisticNonGoalStates, vProbabilistic);
-                storm::utility::vector::setVectorValues(result, markovianNonGoalStates, vMarkovian);
-                return result;
-            }
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        template<typename SparseMarkovAutomatonModelType>
+        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(pathFormula.getLeftSubformula().isTrueFormula(), storm::exceptions::NotImplementedException, "Only bounded properties of the form 'true U[t1, t2] phi' are currently supported.");
-            STORM_LOG_THROW(model.isClosed(), storm::exceptions::InvalidArgumentException, "Unable to compute time-bounded reachability probabilities in non-closed Markov automaton.");
+            STORM_LOG_THROW(this->getModel().isClosed(), storm::exceptions::InvalidArgumentException, "Unable to compute time-bounded reachability probabilities in non-closed Markov automaton.");
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
+            return storm::modelchecker::helper::SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), this->getModel().getMarkovianStates(), rightResult.getTruthValuesVector(), pathFormula.getIntervalBounds());
             std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rightResult.getTruthValuesVector(), pathFormula.getIntervalBounds())));
             return result;
         }
-        
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            return storm::modelchecker::SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(minimize, model.getTransitionMatrix(), model.getBackwardTransitions(), phiStates, psiStates, *MinMaxLinearEquationSolverFactory, qualitative);
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+                
+        template<typename SparseMarkovAutomatonModelType>
+        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
@@ -210,366 +48,34 @@ namespace storm {
             ExplicitQualitativeCheckResult& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative)));
         }
-        
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            std::vector<ValueType> totalRewardVector;
-            if (model.hasTransitionRewards()) {
-                totalRewardVector = model.getTransitionMatrix().getPointwiseProductRowSumVector(model.getTransitionRewardMatrix());
-                if (model.hasStateRewards()) {
-                    storm::utility::vector::addVectors(totalRewardVector, model.getStateRewardVector(), totalRewardVector);
-                }
-            } else {
-                totalRewardVector = std::vector<ValueType>(model.getStateRewardVector());
-            }
-            
-            return this->computeExpectedRewards(minimize, psiStates, totalRewardVector);
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+                
+        template<typename SparseMarkovAutomatonModelType>
+        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(model.isClosed(), storm::exceptions::InvalidArgumentException, "Unable to compute reachability rewards in non-closed Markov automaton.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
         }
-        
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            // If there are no goal states, we avoid the computation and directly return zero.
-            if (psiStates.empty()) {
-                return std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
-            }
-            
-            // Likewise, if all bits are set, we can avoid the computation and set.
-            if ((~psiStates).empty()) {
-                return std::vector<ValueType>(model.getNumberOfStates(), storm::utility::one<ValueType>());
-            }
-            
-            // Start by decomposing the Markov automaton into its MECs.
-            storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(model);
-            
-            // Get some data members for convenience.
-            typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = model.getTransitionMatrix();
-            std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
-            
-            // Now start with compute the long-run average for all end components in isolation.
-            std::vector<ValueType> lraValuesForEndComponents;
-            
-            // While doing so, we already gather some information for the following steps.
-            std::vector<uint_fast64_t> stateToMecIndexMap(model.getNumberOfStates());
-            storm::storage::BitVector statesInMecs(model.getNumberOfStates());
-            
-            for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
-                storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
-                
-                // Gather information for later use.
-                for (auto const& stateChoicesPair : mec) {
-                    uint_fast64_t state = stateChoicesPair.first;
-                    
-                    statesInMecs.set(state);
-                    stateToMecIndexMap[state] = currentMecIndex;
-                }
-                
-                // Compute the LRA value for the current MEC.
-                lraValuesForEndComponents.push_back(this->computeLraForMaximalEndComponent(minimize, transitionMatrix, nondeterministicChoiceIndices, model.getMarkovianStates(), model.getExitRates(), psiStates, mec, currentMecIndex));
-            }
-            
-            // For fast transition rewriting, we build some auxiliary data structures.
-            storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
-            uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
-            uint_fast64_t lastStateNotInMecs = 0;
-            uint_fast64_t numberOfStatesNotInMecs = 0;
-            std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
-            statesNotInMecsBeforeIndex.reserve(model.getNumberOfStates());
-            for (auto state : statesNotContainedInAnyMec) {
-                while (lastStateNotInMecs <= state) {
-                    statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
-                    ++lastStateNotInMecs;
-                }
-                ++numberOfStatesNotInMecs;
-            }
-            
-            // Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
-            std::vector<ValueType> b;
-            typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
-            
-            // If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
-            uint_fast64_t currentChoice = 0;
-            for (auto state : statesNotContainedInAnyMec) {
-                sspMatrixBuilder.newRowGroup(currentChoice);
                 
-                for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
-                    std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
-                    b.push_back(storm::utility::zero<ValueType>());
-                    
-                    for (auto element : transitionMatrix.getRow(choice)) {
-                        if (statesNotContainedInAnyMec.get(element.getColumn())) {
-                            // If the target state is not contained in an MEC, we can copy over the entry.
-                            sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
-                        } else {
-                            // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
-                            // so that we are able to write the cumulative probability to the MEC into the matrix.
-                            auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
-                        }
-                    }
-                    
-                    // Now insert all (cumulative) probability values that target an MEC.
-                    for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
-                        if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
-                            sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
-                        }
-                    }
-                }
-            }
-            
-            // Now we are ready to construct the choices for the auxiliary states.
-            for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
-                storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
-                sspMatrixBuilder.newRowGroup(currentChoice);
-                
-                for (auto const& stateChoicesPair : mec) {
-                    uint_fast64_t state = stateChoicesPair.first;
-                    boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
-                    
-                    for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
-                        std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
-                        
-                        // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
-                        if (choicesInMec.find(choice) == choicesInMec.end()) {
-                            b.push_back(storm::utility::zero<ValueType>());
-                            
-                            for (auto element : transitionMatrix.getRow(choice)) {
-                                if (statesNotContainedInAnyMec.get(element.getColumn())) {
-                                    // If the target state is not contained in an MEC, we can copy over the entry.
-                                    sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
-                                } else {
-                                    // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
-                                    // so that we are able to write the cumulative probability to the MEC into the matrix.
-                                    auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
-                                }
-                            }
-                            
-                            // Now insert all (cumulative) probability values that target an MEC.
-                            for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
-                                if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
-                                    sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
-                                }
-                            }
-                            
-                            ++currentChoice;
-                        }
-                    }
-                }
-                
-                // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
-                ++currentChoice;
-                b.push_back(lraValuesForEndComponents[mecIndex]);
-            }
-            
-            // Finalize the matrix and solve the corresponding system of equations.
-            storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
-            
-            std::vector<ValueType> x(numberOfStatesNotInMecs + mecDecomposition.size());
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(sspMatrix);
-            solver->solveEquationSystem(minimize, x, b);
-            
-            // Prepare result vector.
-            std::vector<ValueType> result(model.getNumberOfStates());
-            
-            // Set the values for states not contained in MECs.
-            storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, x);
-            
-            // Set the values for all states in MECs.
-            for (auto state : statesInMecs) {
-				result[state] = x[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
-            }
-            
-            return result;
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        template<typename SparseMarkovAutomatonModelType>
+        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(model.isClosed(), storm::exceptions::InvalidArgumentException, "Unable to compute long-run average in non-closed Markov automaton.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
         }
-        
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeExpectedTimesHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            std::vector<ValueType> rewardValues(model.getNumberOfStates(), storm::utility::zero<ValueType>());
-            storm::utility::vector::setVectorValues(rewardValues, model.getMarkovianStates(), storm::utility::one<ValueType>());
-            return this->computeExpectedRewards(minimize, psiStates, rewardValues);
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::computeExpectedTimes(storm::logic::EventuallyFormula const& eventuallyFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+                
+        template<typename SparseMarkovAutomatonModelType>
+        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<SparseMarkovAutomatonModelType>::computeExpectedTimes(storm::logic::EventuallyFormula const& eventuallyFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(model.isClosed(), storm::exceptions::InvalidArgumentException, "Unable to compute expected times in non-closed Markov automaton.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(eventuallyFormula.getSubformula());
             ExplicitQualitativeCheckResult& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeExpectedTimesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
         }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::checkBooleanLiteralFormula(storm::logic::BooleanLiteralFormula const& stateFormula) {
-            if (stateFormula.isTrueFormula()) {
-                return std::unique_ptr<CheckResult>(new ExplicitQualitativeCheckResult(storm::storage::BitVector(model.getNumberOfStates(), true)));
-            } else {
-                return std::unique_ptr<CheckResult>(new ExplicitQualitativeCheckResult(storm::storage::BitVector(model.getNumberOfStates())));
-            }
-        }
-        
-        template<typename ValueType>
-        std::unique_ptr<CheckResult> SparseMarkovAutomatonCslModelChecker<ValueType>::checkAtomicLabelFormula(storm::logic::AtomicLabelFormula const& stateFormula) {
-            STORM_LOG_THROW(model.hasLabel(stateFormula.getLabel()), storm::exceptions::InvalidPropertyException, "The property refers to unknown label '" << stateFormula.getLabel() << "'.");
-            return std::unique_ptr<CheckResult>(new ExplicitQualitativeCheckResult(model.getStates(stateFormula.getLabel())));
-        }
                 
-        template<typename ValueType>
-        ValueType SparseMarkovAutomatonCslModelChecker<ValueType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::BitVector const& markovianStates, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec, uint_fast64_t mecIndex) {
-            std::unique_ptr<storm::utility::solver::LpSolverFactory> lpSolverFactory(new storm::utility::solver::LpSolverFactory());
-            std::unique_ptr<storm::solver::LpSolver> solver = lpSolverFactory->create("LRA for MEC");
-            solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
-            
-            // First, we need to create the variables for the problem.
-            std::map<uint_fast64_t, storm::expressions::Variable> stateToVariableMap;
-            for (auto const& stateChoicesPair : mec) {
-                std::string variableName = "x" + std::to_string(stateChoicesPair.first);
-                stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName);
-            }
-            storm::expressions::Variable k = solver->addUnboundedContinuousVariable("k", 1);
-            solver->update();
-            
-            // Now we encode the problem as constraints.
-            for (auto const& stateChoicesPair : mec) {
-                uint_fast64_t state = stateChoicesPair.first;
-                
-                // Now, based on the type of the state, create a suitable constraint.
-                if (markovianStates.get(state)) {
-                    storm::expressions::Expression constraint = stateToVariableMap.at(state);
-                    
-                    for (auto element : transitionMatrix.getRow(nondeterministicChoiceIndices[state])) {
-                        constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
-                    }
-                    
-                    constraint = constraint + solver->getConstant(storm::utility::one<ValueType>() / exitRates[state]) * k;
-                    storm::expressions::Expression rightHandSide = psiStates.get(state) ? solver->getConstant(storm::utility::one<ValueType>() / exitRates[state]) : solver->getConstant(0);
-                    if (minimize) {
-                        constraint = constraint <= rightHandSide;
-                    } else {
-                        constraint = constraint >= rightHandSide;
-                    }
-                    solver->addConstraint("state" + std::to_string(state), constraint);
-                } else {
-                    // For probabilistic states, we want to add the constraint x_s <= sum P(s, a, s') * x_s' where a is the current action
-                    // and the sum ranges over all states s'.
-                    for (auto choice : stateChoicesPair.second) {
-                        storm::expressions::Expression constraint = stateToVariableMap.at(state);
-                        
-                        for (auto element : transitionMatrix.getRow(choice)) {
-                            constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
-                        }
-                        
-                        storm::expressions::Expression rightHandSide = solver->getConstant(storm::utility::zero<ValueType>());
-                        if (minimize) {
-                            constraint = constraint <= rightHandSide;
-                        } else {
-                            constraint = constraint >= rightHandSide;
-                        }
-                        solver->addConstraint("state" + std::to_string(state), constraint);
-                    }
-                }
-            }
-            
-            solver->optimize();
-            return solver->getContinuousValue(k);
-        }
-        
-        template<typename ValueType>
-        std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeExpectedRewards(bool minimize, storm::storage::BitVector const& goalStates, std::vector<ValueType> const& stateRewards) const {
-            // First, we need to check which states have infinite expected time (by definition).
-            storm::storage::BitVector infinityStates;
-            if (minimize) {
-                // If we need to compute the minimum expected times, we have to set the values of those states to infinity that, under all schedulers,
-                // reach a bottom SCC without a goal state.
-                
-                // So we start by computing all bottom SCCs without goal states.
-                storm::storage::StronglyConnectedComponentDecomposition<double> sccDecomposition(model, ~goalStates, true, true);
-                
-                // Now form the union of all these SCCs.
-                storm::storage::BitVector unionOfNonGoalBSccs(model.getNumberOfStates());
-                for (auto const& scc : sccDecomposition) {
-                    for (auto state : scc) {
-                        unionOfNonGoalBSccs.set(state);
-                    }
-                }
-                
-                // Finally, if this union is non-empty, compute the states such that all schedulers reach some state of the union.
-                if (!unionOfNonGoalBSccs.empty()) {
-                    infinityStates = storm::utility::graph::performProbGreater0A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), storm::storage::BitVector(model.getNumberOfStates(), true), unionOfNonGoalBSccs);
-                } else {
-                    // Otherwise, we have no infinity states.
-                    infinityStates = storm::storage::BitVector(model.getNumberOfStates());
-                }
-            } else {
-                // If we maximize the property, the expected time of a state is infinite, if an end-component without any goal state is reachable.
-                
-                // So we start by computing all MECs that have no goal state.
-                storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(model, ~goalStates);
-                
-                // Now we form the union of all states in these end components.
-                storm::storage::BitVector unionOfNonGoalMaximalEndComponents(model.getNumberOfStates());
-                for (auto const& mec : mecDecomposition) {
-                    for (auto const& stateActionPair : mec) {
-                        unionOfNonGoalMaximalEndComponents.set(stateActionPair.first);
-                    }
-                }
-                
-                if (!unionOfNonGoalMaximalEndComponents.empty()) {
-                    // Now we need to check for which states there exists a scheduler that reaches one of the previously computed states.
-                    infinityStates = storm::utility::graph::performProbGreater0E(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), storm::storage::BitVector(model.getNumberOfStates(), true), unionOfNonGoalMaximalEndComponents);
-                } else {
-                    // Otherwise, we have no infinity states.
-                    infinityStates = storm::storage::BitVector(model.getNumberOfStates());
-                }
-            }
-            
-            // Now we identify the states for which values need to be computed.
-            storm::storage::BitVector maybeStates = ~(goalStates | infinityStates);
-            
-            // Then, we can eliminate the rows and columns for all states whose values are already known to be 0.
-            std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
-            storm::storage::SparseMatrix<ValueType> submatrix = model.getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates);
-            
-            // Now prepare the expected reward values for all states so they can be used as the right-hand side of the equation system.
-            std::vector<ValueType> rewardValues(stateRewards);
-            for (auto state : model.getMarkovianStates()) {
-                rewardValues[state] = rewardValues[state] / model.getExitRates()[state];
-            }
-            
-            // Finally, prepare the actual right-hand side.
-            std::vector<ValueType> b(submatrix.getRowCount());
-            storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, model.getNondeterministicChoiceIndices(), rewardValues);
-            
-            // Solve the corresponding system of equations.
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(submatrix);
-            solver->solveEquationSystem(minimize, x, b);
-            
-            // Create resulting vector.
-            std::vector<ValueType> result(model.getNumberOfStates());
-            
-            // Set values of resulting vector according to previous result and return the result.
-            storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
-            storm::utility::vector::setVectorValues(result, goalStates, storm::utility::zero<ValueType>());
-            storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
-            
-            return result;
-        }
-        
         template class SparseMarkovAutomatonCslModelChecker<storm::models::sparse::MarkovAutomaton<double>>;
     }
 }
\ No newline at end of file
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
index cd51608eb..26e47bd92 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
@@ -1,73 +1,35 @@
 #ifndef STORM_MODELCHECKER_CSL_SPARSEMARKOVAUTOMATONCSLMODELCHECKER_H_
 #define STORM_MODELCHECKER_CSL_SPARSEMARKOVAUTOMATONCSLMODELCHECKER_H_
 
-#include "src/modelchecker/AbstractModelChecker.h"
-#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
+#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
+
 #include "src/models/sparse/MarkovAutomaton.h"
-#include "src/storage/MaximalEndComponentDecomposition.h"
-#include "src/solver/MinMaxLinearEquationSolver.h"
+
 #include "src/utility/solver.h"
 
 namespace storm {
     namespace modelchecker {
         
-        template<typename ValueType>
-        class SparseMarkovAutomatonCslModelChecker : public AbstractModelChecker {
+        template<typename SparseMarkovAutomatonModelType>
+        class SparseMarkovAutomatonCslModelChecker : public SparsePropositionalModelChecker<SparseMarkovAutomatonModelType>  {
         public:
-            explicit SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& MinMaxLinearEquationSolver);
-            explicit SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model);
+            typedef typename SparseMarkovAutomatonModelType::ValueType ValueType;
+            typedef typename SparseMarkovAutomatonModelType::RewardModelType RewardModelType;
+            
+            explicit SparseMarkovAutomatonCslModelChecker(SparseMarkovAutomatonModelType const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& minMaxLinearEquationSolver);
+            explicit SparseMarkovAutomatonCslModelChecker(SparseMarkovAutomatonModelType const& model);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(storm::logic::Formula const& formula) const override;
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeExpectedTimes(storm::logic::EventuallyFormula const& eventuallyFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>());
-            virtual std::unique_ptr<CheckResult> checkBooleanLiteralFormula(storm::logic::BooleanLiteralFormula const& stateFormula) override;
-            virtual std::unique_ptr<CheckResult> checkAtomicLabelFormula(storm::logic::AtomicLabelFormula const& stateFormula) override;
             
         private:
-            // The methods that perform the actual checking.
-            std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair) const;
-            static void computeBoundedReachabilityProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory);
-            std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            std::vector<ValueType> computeReachabilityRewardsHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            std::vector<ValueType> computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            std::vector<ValueType> computeExpectedTimesHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
-
-            /*!
-             * Computes the long-run average value for the given maximal end component of a Markov automaton.
-             *
-             * @param minimize Sets whether the long-run average is to be minimized or maximized.
-             * @param transitionMatrix The transition matrix of the underlying Markov automaton.
-             * @param nondeterministicChoiceIndices A vector indicating at which row the choice of a given state begins.
-             * @param markovianStates A bit vector storing all markovian states.
-             * @param exitRates A vector with exit rates for all states. Exit rates of probabilistic states are assumed to be zero.
-             * @param goalStates A bit vector indicating which states are to be considered as goal states.
-             * @param mec The maximal end component to consider for computing the long-run average.
-             * @param mecIndex The index of the MEC.
-             * @return The long-run average of being in a goal state for the given MEC.
-             */
-            static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::BitVector const& markovianStates, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec, uint_fast64_t mecIndex = 0);
-            
-            /*!
-             * Computes the expected reward that is gained from each state before entering any of the goal states.
-             *
-             * @param minimize Indicates whether minimal or maximal rewards are to be computed.
-             * @param goalStates The goal states that define until which point rewards are gained.
-             * @param stateRewards A vector that defines the reward gained in each state. For probabilistic states, this is an instantaneous reward
-             * that is fully gained and for Markovian states the actually gained reward is dependent on the expected time to stay in the
-             * state, i.e. it is gouverned by the exit rate of the state.
-             * @return A vector that contains the expected reward for each state of the model.
-             */
-            std::vector<ValueType> computeExpectedRewards(bool minimize, storm::storage::BitVector const& goalStates, std::vector<ValueType> const& stateRewards) const;
-            
-            // The model this model checker is supposed to analyze.
-            storm::models::sparse::MarkovAutomaton<ValueType> const& model;
-            
             // An object that is used for retrieving solvers for systems of linear equations that are the result of nondeterministic choices.
-            std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> MinMaxLinearEquationSolverFactory;
+            std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> minMaxLinearEquationSolverFactory;
         };
     }
 }
diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp b/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
index e69de29bb..ba43b4060 100644
--- a/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
+++ b/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
@@ -0,0 +1,325 @@
+#include "src/modelchecker/csl/helper/HybridCtmcCslHelper.h"
+
+#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h"
+#include "src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h"
+
+#include "src/storage/dd/CuddAdd.h"
+#include "src/storage/dd/CuddBdd.h"
+#include "src/storage/dd/CuddOdd.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/graph.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                
+                boost::optional<storm::dd::Add<DdType>> modifiedStateRewardVector;
+                if (stateRewardVector) {
+                    modifiedStateRewardVector = stateRewardVector.get() / exitRateVector;
+                }
+                
+                return HybridDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(model, computeProbabilityMatrix(model, rateMatrix, exitRateVector), modifiedStateRewardVector, transitionRewardMatrix, targetStates, qualitative, linearEquationSolverFactory);
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeNextProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& nextStates) {
+                return HybridDtmcPrctlHelper<DdType, ValueType>::computeNextProbabilities(model, computeProbabilityMatrix(model, rateMatrix, exitRateVector), nextStates);
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                return HybridDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(model, computeProbabilityMatrix(model, rateMatrix, exitRateVector), phiStates, psiStates, qualitative, linearEquationSolverFactory);
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                
+                // If the time bounds are [0, inf], we rather call untimed reachability.
+                storm::utility::ConstantsComparator<ValueType> comparator;
+                if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
+                    return computeUntilProbabilities(model, rateMatrix, exitRateVector, phiStates, psiStates, qualitative, linearEquationSolverFactory);
+                }
+                
+                // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0
+                // or t' != inf.
+                
+                // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the
+                // further computations.
+                storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, rateMatrix.notZero(), phiStates, psiStates);
+                STORM_LOG_INFO("Found " << statesWithProbabilityGreater0.getNonZeroCount() << " states with probability greater 0.");
+                storm::dd::Bdd<DdType> statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 && !psiStates;
+                STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNonZeroCount() << " 'maybe' states.");
+                
+                if (!statesWithProbabilityGreater0NonPsi.isZero()) {
+                    if (comparator.isZero(upperBound)) {
+                        // In this case, the interval is of the form [0, 0].
+                        return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+                    } else {
+                        if (comparator.isZero(lowerBound)) {
+                            // In this case, the interval is of the form [0, t].
+                            // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier.
+                            
+                            // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
+                            ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRateVector).getMax();
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                            
+                            // Compute the uniformized matrix.
+                            storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, statesWithProbabilityGreater0NonPsi, uniformizationRate);
+                            
+                            // Compute the vector that is to be added as a compensation for removing the absorbing states.
+                            storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * rateMatrix * psiStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(model.getColumnVariables()) / model.getManager().getConstant(uniformizationRate);
+                            
+                            // Create an ODD for the translation to an explicit representation.
+                            storm::dd::Odd<DdType> odd(statesWithProbabilityGreater0NonPsi);
+                            
+                            // Convert the symbolic parts to their explicit representation.
+                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                            std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
+                            
+                            // Finally compute the transient probabilities.
+                            std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
+                            std::vector<ValueType> subresult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, linearEquationSolverFactory);
+                            
+                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(),
+                                                                                                          (psiStates || !statesWithProbabilityGreater0) && model.getReachableStates(),
+                                                                                                          psiStates.toAdd(), statesWithProbabilityGreater0NonPsi, odd, subresult));
+                        } else if (comparator.isInfinity(upperBound)) {
+                            // In this case, the interval is of the form [t, inf] with t != 0.
+                            
+                            // Start by computing the (unbounded) reachability probabilities of reaching psi states while
+                            // staying in phi states.
+                            std::unique_ptr<CheckResult> unboundedResult = computeUntilProbabilities(model, rateMatrix, exitRateVector, phiStates, psiStates, qualitative, linearEquationSolverFactory);
+                            
+                            // Compute the set of relevant states.
+                            storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
+                            
+                            // Filter the unbounded result such that it only contains values for the relevant states.
+                            unboundedResult->filter(SymbolicQualitativeCheckResult<DdType>(model.getReachableStates(), relevantStates));
+                            
+                            // Build an ODD for the relevant states.
+                            storm::dd::Odd<DdType> odd(relevantStates);
+                            
+                            std::vector<ValueType> result;
+                            if (unboundedResult->isHybridQuantitativeCheckResult()) {
+                                std::unique_ptr<CheckResult> explicitUnboundedResult = unboundedResult->asHybridQuantitativeCheckResult<DdType>().toExplicitQuantitativeCheckResult();
+                                result = std::move(explicitUnboundedResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
+                            } else {
+                                STORM_LOG_THROW(unboundedResult->isSymbolicQuantitativeCheckResult(), storm::exceptions::InvalidStateException, "Expected check result of different type.");
+                                result = unboundedResult->asSymbolicQuantitativeCheckResult<DdType>().getValueVector().template toVector<ValueType>(odd);
+                            }
+                            
+                            // Determine the uniformization rate for the transient probability computation.
+                            ValueType uniformizationRate = 1.02 * (relevantStates.toAdd() * exitRateVector).getMax();
+                            
+                            // Compute the uniformized matrix.
+                            storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, relevantStates, uniformizationRate);
+                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                            
+                            // Compute the transient probabilities.
+                            result = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, linearEquationSolverFactory);
+                            
+                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), !relevantStates && model.getReachableStates(), model.getManager().getAddZero(), relevantStates, odd, result));
+                        } else {
+                            // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
+                            
+                            if (lowerBound != upperBound) {
+                                // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'.
+                                
+                                // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
+                                ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRateVector).getMax();
+                                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                                
+                                // Compute the (first) uniformized matrix.
+                                storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, statesWithProbabilityGreater0NonPsi, uniformizationRate);
+                                
+                                // Create the one-step vector.
+                                storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * rateMatrix * psiStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(model.getColumnVariables()) / model.getManager().getConstant(uniformizationRate);
+                                
+                                // Build an ODD for the relevant states and translate the symbolic parts to their explicit representation.
+                                storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0NonPsi);
+                                storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                                std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
+                                
+                                // Compute the transient probabilities.
+                                std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
+                                std::vector<ValueType> subResult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, linearEquationSolverFactory);
+                                
+                                // Transform the explicit result to a hybrid check result, so we can easily convert it to
+                                // a symbolic qualitative format.
+                                HybridQuantitativeCheckResult<DdType> hybridResult(model.getReachableStates(), psiStates || (!statesWithProbabilityGreater0 && model.getReachableStates()),
+                                                                                   psiStates.toAdd(), statesWithProbabilityGreater0NonPsi, odd, subResult);
+                                
+                                // Compute the set of relevant states.
+                                storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
+                                
+                                // Filter the unbounded result such that it only contains values for the relevant states.
+                                hybridResult.filter(SymbolicQualitativeCheckResult<DdType>(model.getReachableStates(), relevantStates));
+                                
+                                // Build an ODD for the relevant states.
+                                odd = storm::dd::Odd<DdType>(relevantStates);
+                                
+                                std::unique_ptr<CheckResult> explicitResult = hybridResult.toExplicitQuantitativeCheckResult();
+                                std::vector<ValueType> newSubresult = std::move(explicitResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
+                                
+                                // Then compute the transient probabilities of being in such a state after t time units. For this,
+                                // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                                uniformizationRate =  1.02 * (relevantStates.toAdd() * exitRateVector).getMax();
+                                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                                
+                                // If the lower and upper bounds coincide, we have only determined the relevant states at this
+                                // point, but we still need to construct the starting vector.
+                                if (lowerBound == upperBound) {
+                                    odd = storm::dd::Odd<DdType>(relevantStates);
+                                    newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
+                                }
+                                
+                                // Finally, we compute the second set of transient probabilities.
+                                uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, relevantStates, uniformizationRate);
+                                explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                                
+                                newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, linearEquationSolverFactory);
+                                
+                                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), !relevantStates && model.getReachableStates(), model.getManager().getAddZero(), relevantStates, odd, newSubresult));
+                            } else {
+                                // In this case, the interval is of the form [t, t] with t != 0, t != inf.
+                                
+                                // Build an ODD for the relevant states.
+                                storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0);
+                                
+                                std::vector<ValueType> newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
+                                
+                                // Then compute the transient probabilities of being in such a state after t time units. For this,
+                                // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                                ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0.toAdd() * exitRateVector).getMax();
+                                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                                
+                                // Finally, we compute the second set of transient probabilities.
+                                storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, statesWithProbabilityGreater0, uniformizationRate);
+                                storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                                
+                                newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, linearEquationSolverFactory);
+                                
+                                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), !statesWithProbabilityGreater0 && model.getReachableStates(), model.getManager().getAddZero(), statesWithProbabilityGreater0, odd, newSubresult));
+                            }
+                        }
+                    }
+                } else {
+                    return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+                }
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                
+                // Only compute the result if the model has a state-based reward model.
+                STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Create ODD for the translation.
+                storm::dd::Odd<DdType> odd(model.getReachableStates());
+                
+                // Initialize result to state rewards of the model.
+                std::vector<ValueType> result = model.getStateRewardVector().template toVector<ValueType>(odd);
+                
+                // If the time-bound is not zero, we need to perform a transient analysis.
+                if (timeBound > 0) {
+                    ValueType uniformizationRate = 1.02 * exitRateVector.getMax();
+                    STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                    
+                    storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector, model.getReachableStates(), uniformizationRate);
+                    
+                    storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                    result = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, linearEquationSolverFactory);
+                }
+                
+                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, result));
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has a state-based reward model.
+                STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // If the time bound is zero, the result is the constant zero vector.
+                if (timeBound == 0) {
+                    return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getAddZero()));
+                }
+                
+                // Otherwise, we need to perform some computations.
+                
+                // Start with the uniformization.
+                ValueType uniformizationRate = 1.02 * exitRateVector.getMax();
+                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                
+                // Create ODD for the translation.
+                storm::dd::Odd<DdType> odd(model.getReachableStates());
+                
+                // Compute the uniformized matrix.
+                storm::dd::Add<DdType> uniformizedMatrix = computeUniformizedMatrix(model, rateMatrix, exitRateVector,  model.getReachableStates(), uniformizationRate);
+                storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                
+                // Then compute the state reward vector to use in the computation.
+                storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+                if (model.hasTransitionRewards()) {
+                    totalRewardVector += (rateMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+                }
+                std::vector<ValueType> explicitTotalRewardVector = totalRewardVector.template toVector<ValueType>(odd);
+                
+                // Finally, compute the transient probabilities.
+                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, linearEquationSolverFactory);
+                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), std::move(odd), std::move(result)));
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeLongRunAverage(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                storm::dd::Add<DdType> probabilityMatrix = computeProbabilityMatrix(model, rateMatrix, exitRateVector);
+                
+                // Create ODD for the translation.
+                storm::dd::Odd<DdType> odd(model.getReachableStates());
+                
+                storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd);
+                std::vector<ValueType> explicitExitRateVector = exitRateVector.template toVector<ValueType>(odd);
+                
+                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeLongRunAverage(explicitProbabilityMatrix, psiStates.toVector(odd), &explicitExitRateVector, qualitative, linearEquationSolverFactory);
+                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), std::move(odd), std::move(result)));
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            storm::dd::Add<DdType> HybridCtmcCslHelper<DdType, ValueType>::computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate) {
+                STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << ".");
+                STORM_LOG_DEBUG("Keeping " << maybeStates.getNonZeroCount() << " rows.");
+                
+                // Cut all non-maybe rows/columns from the transition matrix.
+                storm::dd::Add<DdType> uniformizedMatrix = transitionMatrix * maybeStates.toAdd() * maybeStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+                
+                // Now perform the uniformization.
+                uniformizedMatrix = uniformizedMatrix / model.getManager().getConstant(uniformizationRate);
+                storm::dd::Add<DdType> diagonal = model.getRowColumnIdentity() * maybeStates.toAdd();
+                storm::dd::Add<DdType> diagonalOffset = diagonal;
+                diagonalOffset -= diagonal * (exitRateVector / model.getManager().getConstant(uniformizationRate));
+                uniformizedMatrix += diagonalOffset;
+                
+                return uniformizedMatrix;
+            }
+            
+            template<storm::dd::DdType DdType, class ValueType>
+            storm::dd::Add<DdType> HybridCtmcCslHelper<DdType, ValueType>::computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector) {
+                return rateMatrix / exitRateVector;
+            }
+
+            template class HybridCtmcCslHelper<storm::dd::DdType::CUDD, double>;
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.h b/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
index e69de29bb..b768ea103 100644
--- a/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
+++ b/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
@@ -0,0 +1,60 @@
+#ifndef STORM_MODELCHECKER_HYBRID_CTMC_CSL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_HYBRID_CTMC_CSL_MODELCHECKER_HELPER_H_
+
+#include <memory>
+
+#include "src/models/symbolic/Ctmc.h"
+
+#include "src/modelchecker/results/CheckResult.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            class HybridCtmcCslHelper {
+            public:
+                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeCumulativeRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeReachabilityRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeLongRunAverage(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeNextProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& nextStates);
+
+                /*!
+                 * Converts the given rate-matrix into a time-abstract probability matrix.
+                 *
+                 * @param model The symbolic model.
+                 * @param rateMatrix The rate matrix.
+                 * @param exitRateVector The exit rate vector of the model.
+                 * @return The probability matrix.
+                 */
+                static storm::dd::Add<DdType> computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector);
+                
+                /*!
+                 * Computes the matrix representing the transitions of the uniformized CTMC.
+                 *
+                 * @param model The symbolic model.
+                 * @param transitionMatrix The matrix to uniformize.
+                 * @param exitRateVector The exit rate vector.
+                 * @param maybeStates The states that need to be considered.
+                 * @param uniformizationRate The rate to be used for uniformization.
+                 * @return The uniformized matrix.
+                 */
+                static storm::dd::Add<DdType> computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_HYBRID_CTMC_CSL_MODELCHECKER_HELPER_H_ */
\ No newline at end of file
diff --git a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
new file mode 100644
index 000000000..089a39c8a
--- /dev/null
+++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -0,0 +1,494 @@
+#include "src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h"
+
+#include "src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h"
+
+#include "src/storage/StronglyConnectedComponentDecomposition.h"
+#include "src/storage/MaximalEndComponentDecomposition.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/vector.h"
+#include "src/utility/graph.h"
+
+#include "src/utility/numerical.h"
+
+#include "src/solver/LpSolver.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            template<typename ValueType, typename RewardModelType>
+            void SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeBoundedReachabilityProbabilities(bool min, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // Start by computing four sparse matrices:
+                // * a matrix aMarkovian with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states.
+                // * a matrix aMarkovianToProbabilistic with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states.
+                // * a matrix aProbabilistic with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states.
+                // * a matrix aProbabilisticToMarkovian with all (non-discretized) transitions from probabilistic non-goal states to all Markovian non-goal states.
+                typename storm::storage::SparseMatrix<ValueType> aMarkovian = transitionMatrix.getSubmatrix(true, markovianNonGoalStates, markovianNonGoalStates, true);
+                typename storm::storage::SparseMatrix<ValueType> aMarkovianToProbabilistic = transitionMatrix.getSubmatrix(true, markovianNonGoalStates, probabilisticNonGoalStates);
+                typename storm::storage::SparseMatrix<ValueType> aProbabilistic = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, probabilisticNonGoalStates);
+                typename storm::storage::SparseMatrix<ValueType> aProbabilisticToMarkovian = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, markovianNonGoalStates);
+                
+                // The matrices with transitions from Markovian states need to be digitized.
+                // Digitize aMarkovian. Based on whether the transition is a self-loop or not, we apply the two digitization rules.
+                uint_fast64_t rowIndex = 0;
+                for (auto state : markovianNonGoalStates) {
+                    for (auto& element : aMarkovian.getRow(rowIndex)) {
+                        ValueType eTerm = std::exp(-exitRates[state] * delta);
+                        if (element.getColumn() == rowIndex) {
+                            element.setValue((storm::utility::one<ValueType>() - eTerm) * element.getValue() + eTerm);
+                        } else {
+                            element.setValue((storm::utility::one<ValueType>() - eTerm) * element.getValue());
+                        }
+                    }
+                    ++rowIndex;
+                }
+                
+                // Digitize aMarkovianToProbabilistic. As there are no self-loops in this case, we only need to apply the digitization formula for regular successors.
+                rowIndex = 0;
+                for (auto state : markovianNonGoalStates) {
+                    for (auto& element : aMarkovianToProbabilistic.getRow(rowIndex)) {
+                        element.setValue((1 - std::exp(-exitRates[state] * delta)) * element.getValue());
+                    }
+                    ++rowIndex;
+                }
+                
+                // Initialize the two vectors that hold the variable one-step probabilities to all target states for probabilistic and Markovian (non-goal) states.
+                std::vector<ValueType> bProbabilistic(aProbabilistic.getRowCount());
+                std::vector<ValueType> bMarkovian(markovianNonGoalStates.getNumberOfSetBits());
+                
+                // Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones.
+                std::vector<ValueType> bProbabilisticFixed = transitionMatrix.getConstrainedRowSumVector(probabilisticNonGoalStates, goalStates);
+                std::vector<ValueType> bMarkovianFixed;
+                bMarkovianFixed.reserve(markovianNonGoalStates.getNumberOfSetBits());
+                for (auto state : markovianNonGoalStates) {
+                    bMarkovianFixed.push_back(storm::utility::zero<ValueType>());
+                    
+                    for (auto& element : transitionMatrix.getRowGroup(state)) {
+                        if (goalStates.get(element.getColumn())) {
+                            bMarkovianFixed.back() += (1 - std::exp(-exitRates[state] * delta)) * element.getValue();
+                        }
+                    }
+                }
+                
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic);
+                
+                // Perform the actual value iteration
+                // * loop until the step bound has been reached
+                // * in the loop:
+                // *    perform value iteration using A_PSwG, v_PS and the vector b where b = (A * 1_G)|PS + A_PStoMS * v_MS
+                //      and 1_G being the characteristic vector for all goal states.
+                // *    perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS
+                std::vector<ValueType> markovianNonGoalValuesSwap(markovianNonGoalValues);
+                std::vector<ValueType> multiplicationResultScratchMemory(aProbabilistic.getRowCount());
+                std::vector<ValueType> aProbabilisticScratchMemory(probabilisticNonGoalValues.size());
+                for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
+                    // Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
+                    aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
+                    storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
+                    
+                    // Now perform the inner value iteration for probabilistic states.
+                    solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+                    
+                    // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
+                    aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
+                    storm::utility::vector::addVectors(bMarkovian, bMarkovianFixed, bMarkovian);
+                    aMarkovian.multiplyWithVector(markovianNonGoalValues, markovianNonGoalValuesSwap);
+                    std::swap(markovianNonGoalValues, markovianNonGoalValuesSwap);
+                    storm::utility::vector::addVectors(markovianNonGoalValues, bMarkovian, markovianNonGoalValues);
+                }
+                
+                // After the loop, perform one more step of the value iteration for PS states.
+                aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
+                storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
+                solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
+                
+                // 'Unpack' the bounds to make them more easily accessible.
+                double lowerBound = boundsPair.first;
+                double upperBound = boundsPair.second;
+                
+                // (1) Compute the accuracy we need to achieve the required error bound.
+                ValueType maxExitRate = 0;
+                for (auto value : exitRateVector) {
+                    maxExitRate = std::max(maxExitRate, value);
+                }
+                ValueType delta = (2 * storm::settings::generalSettings().getPrecision()) / (upperBound * maxExitRate * maxExitRate);
+                
+                // (2) Compute the number of steps we need to make for the interval.
+                uint_fast64_t numberOfSteps = static_cast<uint_fast64_t>(std::ceil((upperBound - lowerBound) / delta));
+                STORM_LOG_INFO("Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [" << lowerBound << ", " << upperBound << "]." << std::endl);
+                
+                // (3) Compute the non-goal states and initialize two vectors
+                // * vProbabilistic holds the probability values of probabilistic non-goal states.
+                // * vMarkovian holds the probability values of Markovian non-goal states.
+                storm::storage::BitVector const& markovianNonGoalStates = markovianStates & ~psiStates;
+                storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~psiStates;
+                std::vector<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
+                std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
+                
+                computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRateVector, markovianStates, psiStates, markovianNonGoalStates, probabilisticNonGoalStates, vMarkovian, vProbabilistic, delta, numberOfSteps, minMaxLinearEquationSolverFactory);
+                
+                // (4) If the lower bound of interval was non-zero, we need to take the current values as the starting values for a subsequent value iteration.
+                if (lowerBound != storm::utility::zero<ValueType>()) {
+                    std::vector<ValueType> vAllProbabilistic((~markovianStates).getNumberOfSetBits());
+                    std::vector<ValueType> vAllMarkovian(markovianStates.getNumberOfSetBits());
+                    
+                    // Create the starting value vectors for the next value iteration based on the results of the previous one.
+                    storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, psiStates % ~markovianStates, storm::utility::one<ValueType>());
+                    storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, ~psiStates % ~markovianStates, vProbabilistic);
+                    storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, psiStates % markovianStates, storm::utility::one<ValueType>());
+                    storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, ~psiStates % markovianStates, vMarkovian);
+                    
+                    // Compute the number of steps to reach the target interval.
+                    numberOfSteps = static_cast<uint_fast64_t>(std::ceil(lowerBound / delta));
+                    STORM_LOG_INFO("Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [0, " << lowerBound << "]." << std::endl);
+                    
+                    // Compute the bounded reachability for interval [0, b-a].
+                    computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRateVector, markovianStates, storm::storage::BitVector(numberOfStates), markovianStates, ~markovianStates, vAllMarkovian, vAllProbabilistic, delta, numberOfSteps, minMaxLinearEquationSolverFactory);
+                    
+                    // Create the result vector out of vAllProbabilistic and vAllMarkovian and return it.
+                    std::vector<ValueType> result(numberOfStates, storm::utility::zero<ValueType>());
+                    storm::utility::vector::setVectorValues(result, ~markovianStates, vAllProbabilistic);
+                    storm::utility::vector::setVectorValues(result, markovianStates, vAllMarkovian);
+                    
+                    return result;
+                } else {
+                    // Create the result vector out of 1_G, vProbabilistic and vMarkovian and return it.
+                    std::vector<ValueType> result(numberOfStates);
+                    storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
+                    storm::utility::vector::setVectorValues(result, probabilisticNonGoalStates, vProbabilistic);
+                    storm::utility::vector::setVectorValues(result, markovianNonGoalStates, vMarkovian);
+                    return result;
+                }
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                return storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType, RewardModelType>::computeUntilProbabilities(minimize, transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, minMaxLinearEquationSolverFactory);
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeReachabilityRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                std::vector<ValueType> totalRewardVector = rewardModel.getTotalRewardVector(transitionMatrix.getRowCount(), transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true));
+                return computeExpectedRewards(minimize, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, psiStates, totalRewardVector, minMaxLinearEquationSolverFactory);
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeLongRunAverage(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
+                
+                // If there are no goal states, we avoid the computation and directly return zero.
+                if (psiStates.empty()) {
+                    return std::vector<ValueType>(numberOfStates, storm::utility::zero<ValueType>());
+                }
+                
+                // Likewise, if all bits are set, we can avoid the computation and set.
+                if ((~psiStates).empty()) {
+                    return std::vector<ValueType>(numberOfStates, storm::utility::one<ValueType>());
+                }
+                
+                // Start by decomposing the Markov automaton into its MECs.
+                storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(transitionMatrix, backwardTransitions);
+                
+                // Get some data members for convenience.
+                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = transitionMatrix.getRowGroupIndices();
+                
+                // Now start with compute the long-run average for all end components in isolation.
+                std::vector<ValueType> lraValuesForEndComponents;
+                
+                // While doing so, we already gather some information for the following steps.
+                std::vector<uint_fast64_t> stateToMecIndexMap(numberOfStates);
+                storm::storage::BitVector statesInMecs(numberOfStates);
+                
+                for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
+                    storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
+                    
+                    // Gather information for later use.
+                    for (auto const& stateChoicesPair : mec) {
+                        uint_fast64_t state = stateChoicesPair.first;
+                        
+                        statesInMecs.set(state);
+                        stateToMecIndexMap[state] = currentMecIndex;
+                    }
+                    
+                    // Compute the LRA value for the current MEC.
+                    lraValuesForEndComponents.push_back(computeLraForMaximalEndComponent(minimize, transitionMatrix, exitRateVector, markovianStates, psiStates, mec));
+                }
+                
+                // For fast transition rewriting, we build some auxiliary data structures.
+                storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
+                uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
+                uint_fast64_t lastStateNotInMecs = 0;
+                uint_fast64_t numberOfStatesNotInMecs = 0;
+                std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
+                statesNotInMecsBeforeIndex.reserve(numberOfStates);
+                for (auto state : statesNotContainedInAnyMec) {
+                    while (lastStateNotInMecs <= state) {
+                        statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
+                        ++lastStateNotInMecs;
+                    }
+                    ++numberOfStatesNotInMecs;
+                }
+                
+                // Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
+                std::vector<ValueType> b;
+                typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
+                
+                // If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
+                uint_fast64_t currentChoice = 0;
+                for (auto state : statesNotContainedInAnyMec) {
+                    sspMatrixBuilder.newRowGroup(currentChoice);
+                    
+                    for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
+                        std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+                        b.push_back(storm::utility::zero<ValueType>());
+                        
+                        for (auto element : transitionMatrix.getRow(choice)) {
+                            if (statesNotContainedInAnyMec.get(element.getColumn())) {
+                                // If the target state is not contained in an MEC, we can copy over the entry.
+                                sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+                            } else {
+                                // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+                                // so that we are able to write the cumulative probability to the MEC into the matrix.
+                                auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+                            }
+                        }
+                        
+                        // Now insert all (cumulative) probability values that target an MEC.
+                        for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
+                            if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
+                                sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
+                            }
+                        }
+                    }
+                }
+                
+                // Now we are ready to construct the choices for the auxiliary states.
+                for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
+                    storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
+                    sspMatrixBuilder.newRowGroup(currentChoice);
+                    
+                    for (auto const& stateChoicesPair : mec) {
+                        uint_fast64_t state = stateChoicesPair.first;
+                        boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
+                        
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                            std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+                            
+                            // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
+                            if (choicesInMec.find(choice) == choicesInMec.end()) {
+                                b.push_back(storm::utility::zero<ValueType>());
+                                
+                                for (auto element : transitionMatrix.getRow(choice)) {
+                                    if (statesNotContainedInAnyMec.get(element.getColumn())) {
+                                        // If the target state is not contained in an MEC, we can copy over the entry.
+                                        sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+                                    } else {
+                                        // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+                                        // so that we are able to write the cumulative probability to the MEC into the matrix.
+                                        auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+                                    }
+                                }
+                                
+                                // Now insert all (cumulative) probability values that target an MEC.
+                                for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
+                                    if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
+                                        sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
+                                    }
+                                }
+                                
+                                ++currentChoice;
+                            }
+                        }
+                    }
+                    
+                    // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
+                    ++currentChoice;
+                    b.push_back(lraValuesForEndComponents[mecIndex]);
+                }
+                
+                // Finalize the matrix and solve the corresponding system of equations.
+                storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
+                
+                std::vector<ValueType> x(numberOfStatesNotInMecs + mecDecomposition.size());
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(sspMatrix);
+                solver->solveEquationSystem(minimize, x, b);
+                
+                // Prepare result vector.
+                std::vector<ValueType> result(numberOfStates);
+                
+                // Set the values for states not contained in MECs.
+                storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, x);
+                
+                // Set the values for all states in MECs.
+                for (auto state : statesInMecs) {
+                    result[state] = x[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
+                }
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeExpectedTimes(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
+                std::vector<ValueType> rewardValues(numberOfStates, storm::utility::zero<ValueType>());
+                storm::utility::vector::setVectorValues(rewardValues, markovianStates, storm::utility::one<ValueType>());
+                return computeExpectedRewards(minimize, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, psiStates, rewardValues, minMaxLinearEquationSolverFactory);
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeExpectedRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, std::vector<ValueType> const& stateRewards, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
+                
+                // First, we need to check which states have infinite expected time (by definition).
+                storm::storage::BitVector infinityStates;
+                if (minimize) {
+                    // If we need to compute the minimum expected times, we have to set the values of those states to infinity that, under all schedulers,
+                    // reach a bottom SCC without a goal state.
+                    
+                    // So we start by computing all bottom SCCs without goal states.
+                    storm::storage::StronglyConnectedComponentDecomposition<double> sccDecomposition(transitionMatrix, ~goalStates, true, true);
+                    
+                    // Now form the union of all these SCCs.
+                    storm::storage::BitVector unionOfNonGoalBSccs(numberOfStates);
+                    for (auto const& scc : sccDecomposition) {
+                        for (auto state : scc) {
+                            unionOfNonGoalBSccs.set(state);
+                        }
+                    }
+                    
+                    // Finally, if this union is non-empty, compute the states such that all schedulers reach some state of the union.
+                    if (!unionOfNonGoalBSccs.empty()) {
+                        infinityStates = storm::utility::graph::performProbGreater0A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, storm::storage::BitVector(numberOfStates, true), unionOfNonGoalBSccs);
+                    } else {
+                        // Otherwise, we have no infinity states.
+                        infinityStates = storm::storage::BitVector(numberOfStates);
+                    }
+                } else {
+                    // If we maximize the property, the expected time of a state is infinite, if an end-component without any goal state is reachable.
+                    
+                    // So we start by computing all MECs that have no goal state.
+                    storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(transitionMatrix, backwardTransitions, ~goalStates);
+                    
+                    // Now we form the union of all states in these end components.
+                    storm::storage::BitVector unionOfNonGoalMaximalEndComponents(numberOfStates);
+                    for (auto const& mec : mecDecomposition) {
+                        for (auto const& stateActionPair : mec) {
+                            unionOfNonGoalMaximalEndComponents.set(stateActionPair.first);
+                        }
+                    }
+                    
+                    if (!unionOfNonGoalMaximalEndComponents.empty()) {
+                        // Now we need to check for which states there exists a scheduler that reaches one of the previously computed states.
+                        infinityStates = storm::utility::graph::performProbGreater0E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, storm::storage::BitVector(numberOfStates, true), unionOfNonGoalMaximalEndComponents);
+                    } else {
+                        // Otherwise, we have no infinity states.
+                        infinityStates = storm::storage::BitVector(numberOfStates);
+                    }
+                }
+                
+                // Now we identify the states for which values need to be computed.
+                storm::storage::BitVector maybeStates = ~(goalStates | infinityStates);
+                
+                // Then, we can eliminate the rows and columns for all states whose values are already known to be 0.
+                std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
+                storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates);
+                
+                // Now prepare the expected reward values for all states so they can be used as the right-hand side of the equation system.
+                std::vector<ValueType> rewardValues(stateRewards);
+                for (auto state : markovianStates) {
+                    rewardValues[state] = rewardValues[state] / exitRateVector[state];
+                }
+                
+                // Finally, prepare the actual right-hand side.
+                std::vector<ValueType> b(submatrix.getRowCount());
+                storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, transitionMatrix.getRowGroupIndices(), rewardValues);
+                
+                // Solve the corresponding system of equations.
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(submatrix);
+                solver->solveEquationSystem(minimize, x, b);
+                
+                // Create resulting vector.
+                std::vector<ValueType> result(numberOfStates);
+                
+                // Set values of resulting vector according to previous result and return the result.
+                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
+                storm::utility::vector::setVectorValues(result, goalStates, storm::utility::zero<ValueType>());
+                storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            ValueType SparseMarkovAutomatonCslHelper<ValueType, RewardModelType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec) {
+                std::unique_ptr<storm::utility::solver::LpSolverFactory> lpSolverFactory(new storm::utility::solver::LpSolverFactory());
+                std::unique_ptr<storm::solver::LpSolver> solver = lpSolverFactory->create("LRA for MEC");
+                solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
+                
+                // First, we need to create the variables for the problem.
+                std::map<uint_fast64_t, storm::expressions::Variable> stateToVariableMap;
+                for (auto const& stateChoicesPair : mec) {
+                    std::string variableName = "x" + std::to_string(stateChoicesPair.first);
+                    stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName);
+                }
+                storm::expressions::Variable k = solver->addUnboundedContinuousVariable("k", 1);
+                solver->update();
+                
+                // Now we encode the problem as constraints.
+                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = transitionMatrix.getRowGroupIndices();
+                for (auto const& stateChoicesPair : mec) {
+                    uint_fast64_t state = stateChoicesPair.first;
+                    
+                    // Now, based on the type of the state, create a suitable constraint.
+                    if (markovianStates.get(state)) {
+                        storm::expressions::Expression constraint = stateToVariableMap.at(state);
+                        
+                        for (auto element : transitionMatrix.getRow(nondeterministicChoiceIndices[state])) {
+                            constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
+                        }
+                        
+                        constraint = constraint + solver->getConstant(storm::utility::one<ValueType>() / exitRateVector[state]) * k;
+                        storm::expressions::Expression rightHandSide = goalStates.get(state) ? solver->getConstant(storm::utility::one<ValueType>() / exitRateVector[state]) : solver->getConstant(0);
+                        if (minimize) {
+                            constraint = constraint <= rightHandSide;
+                        } else {
+                            constraint = constraint >= rightHandSide;
+                        }
+                        solver->addConstraint("state" + std::to_string(state), constraint);
+                    } else {
+                        // For probabilistic states, we want to add the constraint x_s <= sum P(s, a, s') * x_s' where a is the current action
+                        // and the sum ranges over all states s'.
+                        for (auto choice : stateChoicesPair.second) {
+                            storm::expressions::Expression constraint = stateToVariableMap.at(state);
+                            
+                            for (auto element : transitionMatrix.getRow(choice)) {
+                                constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
+                            }
+                            
+                            storm::expressions::Expression rightHandSide = solver->getConstant(storm::utility::zero<ValueType>());
+                            if (minimize) {
+                                constraint = constraint <= rightHandSide;
+                            } else {
+                                constraint = constraint >= rightHandSide;
+                            }
+                            solver->addConstraint("state" + std::to_string(state), constraint);
+                        }
+                    }
+                }
+                
+                solver->optimize();
+                return solver->getContinuousValue(k);
+            }
+            
+            template class SparseMarkovAutomatonCslHelper<double>;
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
new file mode 100644
index 000000000..ec70e97ce
--- /dev/null
+++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
@@ -0,0 +1,68 @@
+#ifndef STORM_MODELCHECKER_SPARSE_MARKOVAUTOMATON_CSL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_MARKOVAUTOMATON_CSL_MODELCHECKER_HELPER_H_
+
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/storage/BitVector.h"
+#include "src/storage/MaximalEndComponent.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+            class SparseMarkovAutomatonCslHelper {
+            public:
+                static std::vector<ValueType> computeBoundedUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeReachabilityRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeLongRunAverage(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeExpectedTimes(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+            private:
+                static void computeBoundedReachabilityProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                /*!
+                 * Computes the long-run average value for the given maximal end component of a Markov automaton.
+                 *
+                 * @param minimize Sets whether the long-run average is to be minimized or maximized.
+                 * @param transitionMatrix The transition matrix of the underlying Markov automaton.
+                 * @param markovianStates A bit vector storing all markovian states.
+                 * @param exitRateVector A vector with exit rates for all states. Exit rates of probabilistic states are 
+                 * assumed to be zero.
+                 * @param goalStates A bit vector indicating which states are to be considered as goal states.
+                 * @param mec The maximal end component to consider for computing the long-run average.
+                 * @return The long-run average of being in a goal state for the given MEC.
+                 */
+                static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec);
+                
+                /*!
+                 * Computes the expected reward that is gained from each state before entering any of the goal states.
+                 *
+                 * @param minimize Indicates whether minimal or maximal rewards are to be computed.
+                 * @param transitionMatrix The transition matrix of the underlying Markov automaton.
+                 * @param backwardTransitions The reversed transition relation of the underlying Markov automaton.
+                 * @param exitRateVector A vector with exit rates for all states. Exit rates of probabilistic states are
+                 * assumed to be zero.
+                 * @param markovianStates A bit vector storing all markovian states.
+                 * @param goalStates The goal states that define until which point rewards are gained.
+                 * @param stateRewards A vector that defines the reward gained in each state. For probabilistic states,
+                 * this is an instantaneous reward that is fully gained and for Markovian states the actually gained
+                 * reward is dependent on the expected time to stay in the state, i.e. it is gouverned by the exit rate
+                 * of the state.
+                 * @return A vector that contains the expected reward for each state of the model.
+                 */
+                static std::vector<ValueType> computeExpectedRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, std::vector<ValueType> const& stateRewards, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_MARKOVAUTOMATON_CSL_MODELCHECKER_HELPER_H_ */
diff --git a/src/storage/MaximalEndComponentDecomposition.cpp b/src/storage/MaximalEndComponentDecomposition.cpp
index cc809e03f..7a8c53cda 100644
--- a/src/storage/MaximalEndComponentDecomposition.cpp
+++ b/src/storage/MaximalEndComponentDecomposition.cpp
@@ -22,6 +22,11 @@ namespace storm {
             performMaximalEndComponentDecomposition(transitionMatrix, backwardTransitions, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true));
         }
         
+        template<typename ValueType>
+        MaximalEndComponentDecomposition<ValueType>::MaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& subsystem) {
+            performMaximalEndComponentDecomposition(transitionMatrix, backwardTransitions, subsystem);
+        }
+        
         template<typename ValueType>
         MaximalEndComponentDecomposition<ValueType>::MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model, storm::storage::BitVector const& subsystem) {
             performMaximalEndComponentDecomposition(model.getTransitionMatrix(), model.getBackwardTransitions(), subsystem);
diff --git a/src/storage/MaximalEndComponentDecomposition.h b/src/storage/MaximalEndComponentDecomposition.h
index 672163e52..9708bd65a 100644
--- a/src/storage/MaximalEndComponentDecomposition.h
+++ b/src/storage/MaximalEndComponentDecomposition.h
@@ -33,6 +33,15 @@ namespace storm  {
              * @param backwardTransition The reversed transition relation.
              */
             MaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions);
+
+            /*
+             * Creates an MEC decomposition of the given subsystem of given model (represented by a row-grouped matrix).
+             *
+             * @param transitionMatrix The transition relation of model to decompose into MECs.
+             * @param backwardTransition The reversed transition relation.
+             * @param subsystem The subsystem to decompose.
+             */
+            MaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& subsystem);
             
             /*!
              * Creates an MEC decomposition of the given subsystem in the given model.