merged master in future branch

Former-commit-id: 0970004f69

CMakeLists.txt

@@ -721,7 +721,12 @@ if (MSVC)
 	target_link_libraries(storm "Dbghelp.lib")   
 endif(MSVC)
 
-# Print Cotire Usage Status
+
+#############################################################
+##
+##	Cotire
+##
+#############################################################
 message (STATUS "StoRM - Using Cotire: ${STORM_USE_COTIRE}")
 
 if (STORM_USE_COTIRE)
@@ -735,12 +740,23 @@ if (STORM_USE_COTIRE)
 	#cotire(storm-performance-tests)
 endif()
 
+#############################################################
+##
+##	libc++abi
+##
+#############################################################
 # Link against libc++abi if requested. May be needed to build on Linux systems using clang.
 if (LINK_LIBCXXABI)
 	message (STATUS "StoRM - Linking against libc++abi.")
 	target_link_libraries(storm "c++abi")
 endif(LINK_LIBCXXABI)
 
+
+#############################################################
+##
+##	Doxygen
+##
+#############################################################
 # Add a target to generate API documentation with Doxygen
 if(DOXYGEN_FOUND)
     set(CMAKE_DOXYGEN_OUTPUT_DIR "${CMAKE_CURRENT_BINARY_DIR}/doc")
@@ -751,6 +767,11 @@ if(DOXYGEN_FOUND)
     add_custom_target(doc ${DOXYGEN_EXECUTABLE} "${CMAKE_CURRENT_BINARY_DIR}/Doxyfile" DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/Doxyfile" COMMENT "Generating API documentation with Doxygen" VERBATIM)
 endif(DOXYGEN_FOUND)
 
+#############################################################
+##
+##	memcheck targets
+##
+#############################################################
 add_custom_target(memcheck valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm -v --fix-deadlocks ${PROJECT_SOURCE_DIR}/examples/dtmc/crowds/crowds5_5.tra examples/dtmc/crowds/crowds5_5.lab DEPENDS storm)
 add_custom_target(memcheck-functional-tests valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm-functional-tests -v --fix-deadlocks	DEPENDS storm-functional-tests)
 add_custom_target(memcheck-performance-tests valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm-performance-tests -v --fix-deadlocks DEPENDS storm-performance-tests)

src/modelchecker/csl/HybridCtmcCslModelChecker.cpp

@@ -141,7 +141,7 @@ namespace storm {
                         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);
+                        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);
@@ -354,7 +354,24 @@ namespace storm {
             
             // Finally, compute the transient probabilities.
             std::vector<ValueType> result = SparseCtmcCslModelChecker<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(), odd, result));
+            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)));
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            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 = SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(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)));
         }
         
         // Explicitly instantiate the model checker.

src/modelchecker/csl/HybridCtmcCslModelChecker.h

@@ -23,7 +23,8 @@ namespace storm {
             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> 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;
             

src/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -12,6 +12,8 @@
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
+#include "src/storage/StronglyConnectedComponentDecomposition.h"
+
 #include "src/exceptions/InvalidStateException.h"
 #include "src/exceptions/InvalidPropertyException.h"
 #include "src/exceptions/NotImplementedException.h"
@@ -189,7 +191,7 @@ namespace storm {
                             std::vector<ValueType> newSubresult = std::vector<ValueType>(relevantStates.getNumberOfSetBits());
                             storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult);
                             storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one<ValueType>());
-
+                            
                             // 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 = storm::utility::zero<ValueType>();
@@ -344,7 +346,7 @@ namespace storm {
                 weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
                 storm::utility::vector::applyPointwise(result, values, result, addAndScale);
             }
-
+            
             return result;
         }
         
@@ -406,7 +408,7 @@ namespace storm {
             }
             uniformizationRate *= 1.02;
             STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-
+            
             storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector());
             
             // Compute the total state reward vector.
@@ -436,6 +438,22 @@ namespace storm {
             return result;
         }
         
+        template<class ValueType>
+        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<ValueType>::computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
+            storm::storage::SparseMatrix<ValueType> generatorMatrix(rateMatrix, true);
+            
+            // Place the negative exit rate on the diagonal.
+            for (uint_fast64_t row = 0; row < generatorMatrix.getRowCount(); ++row) {
+                for (auto& entry : generatorMatrix.getRow(row)) {
+                    if (entry.getColumn() == row) {
+                        entry.setValue(-exitRates[row]);
+                    }
+                }
+            }
+            
+            return generatorMatrix;
+        }
+        
         template<class ValueType>
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
@@ -455,6 +473,240 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative)));
         }
         
+        template<class ValueType>
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
+            ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
+            
+            storm::storage::SparseMatrix<ValueType> probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory)));
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // If there are no goal states, we avoid the computation and directly return zero.
+            uint_fast64_t numOfStates = transitionMatrix.getRowCount();
+            if (psiStates.empty()) {
+                return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
+            }
+            
+            // Likewise, if all bits are set, we can avoid the computation.
+            if (psiStates.full()) {
+                return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
+            }
+            
+            // Start by decomposing the DTMC into its BSCCs.
+            storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowCount(), true), false, true);
+            
+            STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs.");
+            
+            // Get some data members for convenience.
+            ValueType one = storm::utility::one<ValueType>();
+            ValueType zero = storm::utility::zero<ValueType>();
+            
+            // Prepare the vector holding the LRA values for each of the BSCCs.
+            std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
+            
+            // First we check which states are in BSCCs.
+            storm::storage::BitVector statesInBsccs(numOfStates);
+            storm::storage::BitVector firstStatesInBsccs(numOfStates);
+            
+            for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+                storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+                
+                // Gather information for later use.
+                bool first = true;
+                for (auto const& state : bscc) {
+                    statesInBsccs.set(state);
+                    if (first) {
+                        firstStatesInBsccs.set(state);
+                    }
+                    first = false;
+                }
+            }
+            storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
+            
+            STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs.");
+            
+            // Prepare a vector holding the index within all states that are in BSCCs for every state.
+            std::vector<uint_fast64_t> indexInStatesInBsccs;
+            
+            // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC.
+            std::vector<uint_fast64_t> stateToBsccIndexMap;
+            
+            if (!statesInBsccs.empty()) {
+                firstStatesInBsccs = firstStatesInBsccs % statesInBsccs;
+                
+                // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs.
+                storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
+                
+                // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a
+                // multiplication from the right by transposing the system.
+                bsccEquationSystem = bsccEquationSystem.transpose(false, true);
+                
+                // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states.
+                uint_fast64_t lastIndex = 0;
+                uint_fast64_t currentNumberOfSetBits = 0;
+                indexInStatesInBsccs.reserve(transitionMatrix.getRowCount());
+                for (auto index : statesInBsccs) {
+                    while (lastIndex <= index) {
+                        indexInStatesInBsccs.push_back(currentNumberOfSetBits);
+                        ++lastIndex;
+                    }
+                    ++currentNumberOfSetBits;
+                }
+                
+                stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits());
+                for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+                    for (auto const& state : bscc) {
+                        stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex;
+                    }
+                }
+                
+                // Now build the final equation system matrix, the initial guess and the right-hand side in one go.
+                std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
+                storm::storage::SparseMatrixBuilder<ValueType> builder;
+                for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
+                    
+                    // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the
+                    // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the
+                    // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal.
+                    if (firstStatesInBsccs.get(row)) {
+                        uint_fast64_t requiredBscc = stateToBsccIndexMap[row];
+                        storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc];
+                        
+                        for (auto const& state : bscc) {
+                            builder.addNextValue(row, indexInStatesInBsccs[state], one);
+                        }
+                        
+                        bsccEquationSystemRightSide[row] = one;
+                        
+                    } else {
+                        // Otherwise, we copy the row, and subtract 1 from the diagonal.
+                        for (auto& entry : bsccEquationSystem.getRow(row)) {
+                            if (entry.getColumn() == row) {
+                                builder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
+                            } else {
+                                builder.addNextValue(row, entry.getColumn(), entry.getValue());
+                            }
+                        }
+                    }
+                    
+                }
+                
+                // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC.
+                std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero);
+                for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
+                    
+                    for (auto const& state : bscc) {
+                        bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one /  bscc.size();
+                    }
+                }
+                
+                bsccEquationSystem = builder.build();
+                
+                {
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem);
+                    solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
+                }
+                
+                // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour.
+                if (exitRateVector != nullptr) {
+                    std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
+                    for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) {
+                        bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]);
+                    }
+                    
+                    for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) {
+                        bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]];
+                    }
+                }
+                // Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
+                for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
+                    
+                    for (auto const& state : bscc) {
+                        if (psiStates.get(state)) {
+                            bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += bsccEquationSystemSolution[indexInStatesInBsccs[state]];
+                        }
+                    }
+                }
+                
+                for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << ".");
+                }
+            } else {
+                for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
+                    
+                    // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values
+                    // directly (based on whether or not the contained state is a psi state).
+                    if (psiStates.get(*bscc.begin())) {
+                        bsccLra[bsccIndex] = 1;
+                    }
+                }
+                
+                for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << ".");
+                }
+            }
+            
+            std::vector<ValueType> rewardSolution;
+            if (!statesNotInBsccs.empty()) {
+                // Calculate LRA for states not in bsccs as expected reachability rewards.
+                // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a
+                // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability
+                // of the BSCC.
+                
+                std::vector<ValueType> rewardRightSide;
+                rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
+                
+                for (auto state : statesNotInBsccs) {
+                    ValueType reward = zero;
+                    for (auto entry : transitionMatrix.getRow(state)) {
+                        if (statesInBsccs.get(entry.getColumn())) {
+                            reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]];
+                        }
+                    }
+                    rewardRightSide.push_back(reward);
+                }
+                
+                storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
+                rewardEquationSystemMatrix.convertToEquationSystem();
+                
+                rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one);
+                
+                {
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix);
+                    solver->solveEquationSystem(rewardSolution, rewardRightSide);
+                }
+            }
+            
+            // Fill the result vector.
+            std::vector<ValueType> result(numOfStates);
+            auto rewardSolutionIter = rewardSolution.begin();
+            
+            for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
+                
+                for (auto const& state : bscc) {
+                    result[state] = bsccLra[bsccIndex];
+                }
+            }
+            for (auto state : statesNotInBsccs) {
+                STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution.");
+                // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th
+                // entry in rewardSolution we can just take the next value from the iterator.
+                result[state] = *rewardSolutionIter;
+                ++rewardSolutionIter;
+            }
+            
+            return result;
+        }
+        
+        
         // Explicitly instantiate the model checker.
         template class SparseCtmcCslModelChecker<double>;
         

src/modelchecker/csl/SparseCtmcCslModelChecker.h

@@ -13,10 +13,14 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         class HybridCtmcCslModelChecker;
         
+        template<typename ValueType>
+        class SparseDtmcPrctlModelChecker;
+        
         template<class ValueType>
         class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
             friend class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, ValueType>;
+            friend class SparseDtmcPrctlModelChecker<ValueType>;
             
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model);
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
@@ -29,6 +33,7 @@ namespace storm {
             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> 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::sparse::Ctmc<ValueType> const& getModel() const override;
@@ -39,7 +44,7 @@ namespace storm {
             static std::vector<ValueType> computeUntilProbabilitiesHelper(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::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             std::vector<ValueType> computeInstantaneousRewardsHelper(double timeBound) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const;
-
+            static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             /*!
              * Computes the matrix representing the transitions of the uniformized CTMC.
              *
@@ -73,8 +78,18 @@ namespace storm {
              *
              * @param rateMatrix The rate matrix.
              * @param exitRates The exit rate vector.
+             * @return The †ransition matrix of the embedded DTMC.
              */
             static storm::storage::SparseMatrix<ValueType> computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates);
+
+            /*!
+             * Converts the given rate-matrix into the generator matrix.
+             *
+             * @param rateMatrix The rate matrix.
+             * @param exitRates The exit rate vector.
+             * @return The generator matrix.
+             */
+            static storm::storage::SparseMatrix<ValueType> computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates);
             
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;

src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp

@@ -1,4 +1,5 @@
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 
 #include "src/storage/dd/CuddOdd.h"
 
@@ -301,6 +302,20 @@ namespace storm {
             return this->template getModelAs<storm::models::symbolic::Dtmc<DdType>>();
         }
         
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridDtmcPrctlModelChecker<DdType, ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            
+            // Create ODD for the translation.
+            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
+            
+            storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = this->getModel().getTransitionMatrix().toMatrix(odd, odd);
+            
+            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), 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)));
+        }
+        
         template class HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double>;
     }
 }

src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h

@@ -26,6 +26,7 @@ namespace storm {
             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> 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> 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> 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::Dtmc<DdType> const& getModel() const override;

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@@ -8,11 +8,11 @@
 #include "src/utility/graph.h"
 #include "src/utility/solver.h"
 
+#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
 #include "src/exceptions/InvalidStateException.h"
-#include "src/storage/StronglyConnectedComponentDecomposition.h"
 
 #include "src/exceptions/InvalidPropertyException.h"
 
@@ -302,232 +302,18 @@ namespace storm {
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative)));
         }
-        
-		template<typename ValueType>
-		std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const {
-			// If there are no goal states, we avoid the computation and directly return zero.
-			auto numOfStates = this->getModel().getNumberOfStates();
-			if (psiStates.empty()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
-			}
-
-			// Likewise, if all bits are set, we can avoid the computation and set.
-			if ((~psiStates).empty()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
-			}
-
-			// Start by decomposing the DTMC into its BSCCs.
-			storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(this->getModel(), false, true);
-
-			// Get some data members for convenience.
-			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
-
-			// First we check which states are in BSCCs. We use this later to speed up reachability analysis
-			storm::storage::BitVector statesInBsccs(numOfStates);
-			
-			std::vector<uint_fast64_t> stateToBsccIndexMap(transitionMatrix.getColumnCount());
-
-			for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-				storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-
-				// Gather information for later use.
-				for (auto const& state : bscc) {
-					statesInBsccs.set(state);
-					stateToBsccIndexMap[state] = currentBsccIndex;
-				}
-			}
-
-			storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
-
-			// calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of the transposed transition matrix for the bsccs
-			storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
-
-			//subtract identity matrix
-			for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
-				for (auto& entry : bsccEquationSystem.getRow(row)) {
-					if (entry.getColumn() == row) {
-						entry.setValue(entry.getValue() - one);
-					}
-				}
-			}
-			//now transpose, this internally removes all explicit zeros from the matrix that where introduced when subtracting the identity matrix
-			bsccEquationSystem = bsccEquationSystem.transpose();
-
-			std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
-			std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one);
-			{
-				auto solver = this->linearEquationSolverFactory->create(bsccEquationSystem);
-				solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
-			}
-
-			//calculate LRA Value for each BSCC from steady state distribution in BSCCs
-			// we have to do some scaling, as the probabilities for each BSCC have to sum up to one, which they don't necessarily do in the solution of the equation system
-			std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
-			std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
-			size_t i = 0;
-			for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
-				if (psiStates.get(*stateIter)) {
-					bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
-				}
-				bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
-			}
-			for (i = 0; i < bsccLra.size(); ++i) {
-				bsccLra[i] /= bsccTotalValue[i];
-			}
-
-			//calculate LRA for states not in bsccs as expected reachability rewards
-			//target states are states in bsccs, transition reward is the lra of the bscc for each transition into a bscc and 0 otherwise
-			//this corresponds to sum of LRAs in BSCC weighted by the reachability probability of the BSCC
-
-			std::vector<ValueType> rewardRightSide;
-			rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
-
-			for (auto state : statesNotInBsccs) {
-				ValueType reward = zero;
-				for (auto entry : transitionMatrix.getRow(state)) {
-					if (statesInBsccs.get(entry.getColumn())) {
-						reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
-					}
-				}
-				rewardRightSide.push_back(reward);
-			}
-
-			storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
-			rewardEquationSystemMatrix.convertToEquationSystem();
-
-			std::vector<ValueType> rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one);
-
-			{
-				auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix);
-				solver->solveEquationSystem(rewardSolution, rewardRightSide);
-			}
-
-			// now fill the result vector
-			std::vector<ValueType> result(numOfStates);
-
-			auto rewardSolutionIter = rewardSolution.begin();
-			for (size_t state = 0; state < numOfStates; ++state) {
-				if (statesInBsccs.get(state)) {
-					//assign the value of the bscc the state is in
-					result[state] = bsccLra[stateToBsccIndexMap[state]];
-				} else {
-					assert(rewardSolutionIter != rewardSolution.end());
-					//take the value from the reward computation
-					//since the n-th state not in any bscc is the n-th entry in rewardSolution we can just take the next value from the iterator
-					result[state] = *rewardSolutionIter;
-					++rewardSolutionIter;
-				}
-			}
-
-			return result;
-
-			//old implementeation
-
-			//now we calculate the long run average for each BSCC in isolation and compute the weighted contribution of the BSCC to the LRA value of all states
-			//for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-			//	storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-
-			//	storm::storage::BitVector statesInThisBscc(numOfStates);
-			//	for (auto const& state : bscc) {
-			//		statesInThisBscc.set(state);
-			//	}
-
-			//	//ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc);
-			//	ValueType lraForThisBscc = bsccLra[currentBsccIndex];
-
-			//	//the LRA value of a BSCC contributes to the LRA value of a state with the probability of reaching that BSCC from that state
-			//	//thus we add Prob(Eventually statesInThisBscc) * lraForThisBscc to the result vector
-
-			//	//the reachability probabilities will be zero in other BSCCs, thus we can set the left operand of the until formula to statesNotInBsccs as an optimization
-			//	std::vector<ValueType> reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false);
-			//	
-			//	storm::utility::vector::scaleVectorInPlace<ValueType>(reachProbThisBscc, lraForThisBscc);
-			//	storm::utility::vector::addVectorsInPlace<ValueType>(result, reachProbThisBscc);
-			//}
-
-			//return result;
-		}
-
-		template<typename ValueType>
-		std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-			std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
-			ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-
-			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative)));
-		}
-
 
+        template<typename ValueType>
+        std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
+            ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
+            
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
+        }
+        
 		template<typename ValueType>
-		ValueType SparseDtmcPrctlModelChecker<ValueType>::computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-			//if size is one we can immediately derive the result
-			if (bscc.size() == 1){
-				if (psiStates.get(*(bscc.begin()))) {
-					return storm::utility::one<ValueType>();
-				} else{
-					return storm::utility::zero<ValueType>();
-				}
-			}
-
-			storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount());
-			subsystem.set(bscc.begin(), bscc.end());
-
-			//we now have to solve ((P')^t - I) * x = 0, where P' is the submatrix of transitionMatrix,
-			// ^t means transose, and I is the identity matrix.
-			
-			storm::storage::SparseMatrix<ValueType> subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem, true);
-			subsystemMatrix = subsystemMatrix.transpose();
-
-			// subtract 1 on the diagonal and at the same time add a row with all ones to enforce that the result of the equation system is unique
-			storm::storage::SparseMatrixBuilder<ValueType> equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount());
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
-			bool foundDiagonalElement = false;
-			for (uint_fast64_t row = 0; row < subsystemMatrix.getRowCount(); ++row) {
-				for (auto& entry : subsystemMatrix.getRowGroup(row)) {
-					if (entry.getColumn() == row) {
-						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
-						foundDiagonalElement = true;
-					} else {
-						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue());
-					}
-				}
-
-				// Throw an exception if a row did not have an element on the diagonal.
-				STORM_LOG_THROW(foundDiagonalElement, storm::exceptions::InvalidOperationException, "Internal Error, no diagonal entry found.");
-			}
-			for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) {
-				equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one);
-			}
-			equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero);
-			subsystemMatrix = equationSystemBuilder.build();
-
-			// create x and b for the equation system. setting the last entry of b to one enforces the sum over the unique solution vector is one
-			std::vector<ValueType> steadyStateDist(subsystemMatrix.getRowCount(), zero);
-			std::vector<ValueType> b(subsystemMatrix.getRowCount(), zero);
-			b[subsystemMatrix.getRowCount() - 1] = one;
-
-
-			{
-				auto solver = linearEquationSolverFactory.create(subsystemMatrix);
-				solver->solveEquationSystem(steadyStateDist, b);
-			}
-
-			//remove the last entry of the vector as it was just there to enforce the unique solution
-			steadyStateDist.pop_back();
-			
-			//calculate the fraction we spend in psi states on the long run
-			std::vector<ValueType> steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero);
-			storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist);
-
-			ValueType result = zero;
-
-			for (auto value : steadyStateForPsiStates) {
-				result += value;
-			}
-
-			return result;
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            return SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory);
 		}
 
         template<typename ValueType>

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@@ -9,12 +9,16 @@
 
 namespace storm {
     namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class HybridDtmcPrctlModelChecker;
+        
         // Forward-declare CTMC model checker so we can make it a friend.
         template<typename ValueType> class SparseCtmcCslModelChecker;
         
         template<class ValueType>
         class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
+            friend class HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, ValueType>;
             friend class SparseCtmcCslModelChecker<ValueType>;
             
             explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model);
@@ -41,9 +45,7 @@ namespace storm {
             std::vector<ValueType> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const;
             static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
-			std::vector<ValueType> computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const;
-
-			static ValueType computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
 
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;

src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp

@@ -0,0 +1,263 @@
+#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.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/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        bool SymbolicDtmcPrctlModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const {
+            return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 and 1 of satisfying the until-formula.
+            std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates);
+            storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
+            
+            // Perform some logging.
+            STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
+            STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact probabilities for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                    // maybe states.
+                    storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
+                    prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
+                    subvector = subvector.sumAbstract(model.getColumnVariables());
+                    
+                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                    // for solving the equation system (i.e. compute (I-A)).
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                    
+                    // Solve the equation system.
+                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                    
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + result));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<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());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<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(), this->computeNextProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector())));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
+            storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+            return result.sumAbstract(model.getColumnVariables());
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeBoundedUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 or 1 of satisfying the until-formula.
+            storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound);
+            storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
+            
+            // If there are maybe states, we need to perform matrix-vector multiplications.
+            if (!maybeStates.isZero()) {
+                // Create the ODD for the translation between symbolic and explicit storage.
+                storm::dd::Odd<DdType> odd(maybeStates);
+                
+                // Create the matrix and the vector for the equation system.
+                storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                
+                // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                // non-maybe states in the matrix.
+                storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                
+                // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                // maybe states.
+                storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
+                
+                // Finally cut away all columns targeting non-maybe states.
+                submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &subvector, stepBound);
+                
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd() + result));
+            } else {
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeCumulativeRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Compute the reward vector to add in each step based on the available reward models.
+            storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+            if (model.hasTransitionRewards()) {
+                totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+            }
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &totalRewardVector, stepBound);
+            
+            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeInstantaneousRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getStateRewardVector(), nullptr, stepBound);
+            
+            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, 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>();
+            return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative) {
+            
+            // Only compute the result if there is at least one reward model.
+            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Determine which states have a reward of infinity by definition.
+            storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates);
+            infinityStates = !infinityStates && model.getReachableStates();
+            storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
+            STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
+            STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact rewards for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 1 to indicate that their reward values
+                // are neither 0 nor infinity.
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the state reward vector to use in the computation.
+                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                    if (transitionRewardMatrix) {
+                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
+                    }
+                    
+                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                    // for solving the equation system (i.e. compute (I-A)).
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                    
+                    // Solve the equation system.
+                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                    
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + result));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::models::symbolic::Dtmc<DdType> const& SymbolicDtmcPrctlModelChecker<DdType, ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::symbolic::Dtmc<DdType>>();
+        }
+        
+        template class SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double>;
+    }
+}

src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h

@@ -0,0 +1,49 @@
+#ifndef STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicCtmcCslModelChecker;
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
+        public:
+            friend class SymbolicCtmcCslModelChecker<DdType, ValueType>;
+            
+            explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model);
+            explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory);
+            
+            // 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> 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> 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> 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> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            
+        protected:
+            storm::models::symbolic::Dtmc<DdType> const& getModel() const override;
+            
+        private:
+            // The methods that perform the actual checking.
+            static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static storm::dd::Add<DdType> computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
+            static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative);
+            
+            // An object that is used for retrieving linear equation solvers.
+            std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>> linearEquationSolverFactory;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_ */

src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp

@@ -36,6 +36,11 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new SymbolicQualitativeCheckResult<Type>(model.getReachableStates(), model.getStates(stateFormula.getLabel())));
         }
         
+        template<storm::dd::DdType Type>
+        std::unique_ptr<CheckResult> SymbolicPropositionalModelChecker<Type>::checkAtomicExpressionFormula(storm::logic::AtomicExpressionFormula const& stateFormula) {
+            return std::unique_ptr<CheckResult>(new SymbolicQualitativeCheckResult<Type>(model.getReachableStates(), model.getStates(stateFormula.getExpression())));
+        }
+        
         template<storm::dd::DdType Type>
         storm::models::symbolic::Model<Type> const& SymbolicPropositionalModelChecker<Type>::getModel() const {
             return model;

src/modelchecker/propositional/SymbolicPropositionalModelChecker.h

@@ -16,7 +16,8 @@ namespace storm {
             virtual bool canHandle(storm::logic::Formula const& formula) const override;
             virtual std::unique_ptr<CheckResult> checkBooleanLiteralFormula(storm::logic::BooleanLiteralFormula const& stateFormula) override;
             virtual std::unique_ptr<CheckResult> checkAtomicLabelFormula(storm::logic::AtomicLabelFormula const& stateFormula) override;
-            
+            virtual std::unique_ptr<CheckResult> checkAtomicExpressionFormula(storm::logic::AtomicExpressionFormula const& stateFormula) override;
+
         protected:
             /*!
              * Retrieves the model associated with this model checker instance.

src/modelchecker/reachability/CollectConstraints.h

@@ -1,63 +0,0 @@
-/** 
- * @file:   CollectConstraints.h
- * @author: Sebastian Junges
- *
- * @since October 8, 2014
- */
-
-#pragma once
-#include "src/models/Dtmc.h"
-
-namespace storm {
-namespace modelchecker {
-	namespace reachability {
-		template<typename ValueType>
-		class CollectConstraints
-		{
-			private:
-				std::unordered_set<carl::Constraint<ValueType>> wellformedConstraintSet;
-				std::unordered_set<carl::Constraint<ValueType>> graphPreservingConstraintSet;
-            storm::utility::ConstantsComparator<ValueType> comparator;
-
-			public:
-				std::unordered_set<carl::Constraint<ValueType>> const&  wellformedConstraints() const {
-					return this->wellformedConstraintSet;
-				}
-				
-				std::unordered_set<carl::Constraint<ValueType>> const&  graphPreservingConstraints() const {
-					return this->graphPreservingConstraintSet;
-				}
-				
-				void process(storm::models::Dtmc<ValueType> const& dtmc) 
-				{
-					for(uint_fast64_t state = 0; state < dtmc.getNumberOfStates(); ++state)
-					{
-						ValueType sum;
-						assert(comparator.isZero(sum));
-						for(auto const& transition : dtmc.getRows(state))
-						{
-							sum += transition.getValue();
-							if(!transition.getValue().isConstant())
-							{
-								wellformedConstraintSet.emplace(transition.getValue() - 1, storm::CompareRelation::LEQ);
-								wellformedConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GEQ);
-								graphPreservingConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GT);
-							}
-						}
-						assert(!comparator.isConstant(sum) || comparator.isOne(sum));
-						if(!sum.isConstant()) {
-							wellformedConstraintSet.emplace(sum - 1, storm::CompareRelation::EQ);
-						}
-
-					}
-				}
-
-				void operator()(storm::models::Dtmc<ValueType> const& dtmc) 
-				{
-					process(dtmc);
-				}
-
-		};	
-	}
-}
-}

src/modelchecker/reachability/DirectEncoding.h

@@ -1,80 +0,0 @@
-/**
- * @file:   DirectEncoding.h
- * @author: Sebastian Junges
- *
- * @since April 8, 2014
- */
-
-#pragma once
-
-#ifdef STORM_HAVE_CARL
-#include <carl/io/WriteTosmt2Stream.h>
-
-namespace storm
-{
-    namespace modelchecker
-    {
-        namespace reachability
-        {
-            class DirectEncoding
-            {
-            public:
-                template<typename T>
-                std::string encodeAsSmt2(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& oneStepProbabilities, std::set<carl::Variable> const& parameters, storm::storage::BitVector const& initialStates, typename T::CoeffType const& threshold, bool lessequal = false) {
-                    
-                    carl::io::WriteTosmt2Stream smt2;
-                    uint_fast64_t numberOfStates = transitionMatrix.getRowCount();
-                    carl::VariablePool& vpool = carl::VariablePool::getInstance();
-                    std::vector<carl::Variable> stateVars;
-                    for (carl::Variable const& p : parameters) {
-                        smt2 << ("parameter_bound_" + vpool.getName(p));
-                        smt2 << carl::io::smt2node::AND;
-                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(p), carl::CompareRelation::GT);
-                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(p) - Polynomial::PolyType(1), carl::CompareRelation::LT);
-                        smt2 << carl::io::smt2node::CLOSENODE;
-                    }
-                    
-                    for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
-                        carl::Variable stateVar = vpool.getFreshVariable("s_" + std::to_string(state));
-                        stateVars.push_back(stateVar);
-                        smt2 << ("state_bound_" + std::to_string(state));
-                        smt2 << carl::io::smt2node::AND;
-                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar), carl::CompareRelation::GT);
-                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar) - Polynomial::PolyType(1), carl::CompareRelation::LT);
-                        smt2 << carl::io::smt2node::CLOSENODE;
-                    }
-                    
-                    smt2.setAutomaticLineBreaks(true);
-                    Polynomial::PolyType initStateReachSum;
-                    for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
-                        T reachpropPol;
-                        for (auto const& transition : transitionMatrix.getRow(state)) {
-//                            reachpropPol += transition.getValue() * stateVars[transition.getColumn()];
-                        }
-                        reachpropPol += oneStepProbabilities[state];
-                        smt2 << ("transition_" + std::to_string(state));
-//                        smt2 << carl::Constraint<Polynomial::PolyType>(reachpropPol - stateVars[state], carl::CompareRelation::EQ);
-                    }
-                    
-                    smt2 << ("reachability");
-                    
-                    carl::CompareRelation thresholdRelation = lessequal ? carl::CompareRelation::LEQ : carl::CompareRelation::GEQ;
-                    smt2 << carl::io::smt2node::OR;
-                    for (uint_fast64_t state : initialStates) {
-                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVars[state]) - threshold, thresholdRelation);
-                    }
-                    smt2 << carl::io::smt2node::CLOSENODE;
-                    
-                    smt2 << carl::io::smt2flag::CHECKSAT;
-                    smt2 << carl::io::smt2flag::MODEL;
-                    smt2 << carl::io::smt2flag::UNSAT_CORE;
-                    std::stringstream strm;
-                    strm << smt2;
-                    return strm.str();
-                }
-            };
-        }
-    }
-}
-
-#endif

src/modelchecker/results/ExplicitQualitativeCheckResult.cpp

@@ -29,13 +29,19 @@ namespace storm {
             // Intentionally left empty.
         }
         
-        void ExplicitQualitativeCheckResult::performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, std::function<bool (bool, bool)> const& function) {
+        void ExplicitQualitativeCheckResult::performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, bool logicalAnd) {
             STORM_LOG_THROW(second.isExplicitQualitativeCheckResult(), storm::exceptions::InvalidOperationException, "Cannot perform logical 'and' on check results of incompatible type.");
             STORM_LOG_THROW(first.isResultForAllStates() == second.isResultForAllStates(), storm::exceptions::InvalidOperationException, "Cannot perform logical 'and' on check results of incompatible type.");
             ExplicitQualitativeCheckResult const& secondCheckResult = static_cast<ExplicitQualitativeCheckResult const&>(second);
             if (first.isResultForAllStates()) {
-                boost::get<vector_type>(first.truthValues) &= boost::get<vector_type>(secondCheckResult.truthValues);
+                if (logicalAnd) {
+                    boost::get<vector_type>(first.truthValues) &= boost::get<vector_type>(secondCheckResult.truthValues);
+                } else {
+                    boost::get<vector_type>(first.truthValues) |= boost::get<vector_type>(secondCheckResult.truthValues);
+                }
             } else {
+                std::function<bool (bool, bool)> function = logicalAnd ? [] (bool a, bool b) { return a && b; } : [] (bool a, bool b) { return a || b; };
+                
                 map_type& map1 = boost::get<map_type>(first.truthValues);
                 map_type const& map2 = boost::get<map_type>(secondCheckResult.truthValues);
                 for (auto& element1 : map1) {
@@ -53,12 +59,12 @@ namespace storm {
         }
         
         QualitativeCheckResult& ExplicitQualitativeCheckResult::operator&=(QualitativeCheckResult const& other) {
-            performLogicalOperation(*this, other, [] (bool a, bool b) { return a && b; });
+            performLogicalOperation(*this, other, true);
             return *this;
         }
         
         QualitativeCheckResult& ExplicitQualitativeCheckResult::operator|=(QualitativeCheckResult const& other) {
-            performLogicalOperation(*this, other, [] (bool a, bool b) { return a || b; });
+            performLogicalOperation(*this, other, false);
             return *this;
         }
         

src/modelchecker/results/ExplicitQualitativeCheckResult.h

@@ -51,7 +51,7 @@ namespace storm {
             virtual void filter(QualitativeCheckResult const& filter) override;
 
         private:
-            static void performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, std::function<bool (bool, bool)> const& function);
+            static void performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, bool logicalAnd);
             
             // The values of the quantitative check result.
             boost::variant<vector_type, map_type> truthValues;

src/models/ModelType.cpp

@@ -8,6 +8,7 @@ namespace storm {
                 case ModelType::Ctmc: os << "CTMC"; break;
                 case ModelType::Mdp: os << "MDP"; break;
                 case ModelType::MarkovAutomaton: os << "Markov Automaton"; break;
+                case ModelType::S2pg: os << "S2PG"; break;
             }
             return os;
         }

src/models/ModelType.h

@@ -7,7 +7,7 @@ namespace storm {
     namespace models {
         // All supported model types.
         enum class ModelType {
-            Dtmc, Ctmc, Mdp, MarkovAutomaton
+            Dtmc, Ctmc, Mdp, MarkovAutomaton, S2pg
         };
         
         std::ostream& operator<<(std::ostream& os, ModelType const& type);

src/models/sparse/Ctmc.cpp

@@ -10,7 +10,7 @@ namespace storm {
             Ctmc<ValueType>::Ctmc(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::models::sparse::StateLabeling const& stateLabeling,
                                   boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                   boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                  boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                  boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : DeterministicModel<ValueType>(storm::models::ModelType::Ctmc, rateMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
                 exitRates = createExitRateVector(this->getTransitionMatrix());
             }
@@ -19,7 +19,7 @@ namespace storm {
             Ctmc<ValueType>::Ctmc(storm::storage::SparseMatrix<ValueType>&& rateMatrix, storm::models::sparse::StateLabeling&& stateLabeling,
                                   boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                   boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                  boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                  boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : DeterministicModel<ValueType>(storm::models::ModelType::Ctmc, std::move(rateMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) {
                 // It is important to refer to the transition matrix here, because the given rate matrix has been move elsewhere.
                 exitRates = createExitRateVector(this->getTransitionMatrix());

src/models/sparse/Ctmc.h

@@ -29,7 +29,7 @@ namespace storm {
                 Ctmc(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::models::sparse::StateLabeling const& stateLabeling,
                      boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                      boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                     boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                     boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -43,7 +43,7 @@ namespace storm {
                 Ctmc(storm::storage::SparseMatrix<ValueType>&& rateMatrix, storm::models::sparse::StateLabeling&& stateLabeling,
                      boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                      boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                     boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                     boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 Ctmc(Ctmc<ValueType> const& ctmc) = default;
                 Ctmc& operator=(Ctmc<ValueType> const& ctmc) = default;

src/models/sparse/DeterministicModel.cpp

@@ -11,7 +11,7 @@ namespace storm {
                                                               storm::models::sparse::StateLabeling const& stateLabeling,
                                                               boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                                               boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                                              boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                                              boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : Model<ValueType>(modelType, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
                 // Intentionally left empty.
             }
@@ -22,7 +22,7 @@ namespace storm {
                                                               storm::models::sparse::StateLabeling&& stateLabeling,
                                                               boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                                               boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                                              boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                                              boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : Model<ValueType>(modelType, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) {
                 // Intentionally left empty.
             }

src/models/sparse/DeterministicModel.h

@@ -29,7 +29,7 @@ namespace storm {
                                    storm::models::sparse::StateLabeling const& stateLabeling,
                                    boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                    boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                   boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                   boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -46,7 +46,7 @@ namespace storm {
                                    storm::models::sparse::StateLabeling&& stateLabeling,
                                    boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                    boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                   boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                   boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 DeterministicModel(DeterministicModel const& other) = default;
                 DeterministicModel& operator=(DeterministicModel const& other) = default;

src/models/sparse/Dtmc.cpp

@@ -12,7 +12,7 @@ namespace storm {
                  storm::models::sparse::StateLabeling const& stateLabeling,
                  boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                  boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                 boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                 boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : DeterministicModel<ValueType>(storm::models::ModelType::Dtmc, probabilityMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
                 STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid.");
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
@@ -22,7 +22,7 @@ namespace storm {
             Dtmc<ValueType>::Dtmc(storm::storage::SparseMatrix<ValueType>&& probabilityMatrix, storm::models::sparse::StateLabeling&& stateLabeling,
                  boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                  boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                 boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                 boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : DeterministicModel<ValueType>(storm::models::ModelType::Dtmc, std::move(probabilityMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) {
                 STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid.");
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
@@ -41,7 +41,7 @@ namespace storm {
                 //					  	  	  	  	  	  storm::models::sparse::StateLabeling(this->getStateLabeling(), subSysStates),
                 //					  	  	  	  	  	  boost::optional<std::vector<ValueType>>(),
                 //					  	  	  	  	  	  boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                //					  	  	  	  	  	  boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                //					  	  	  	  	  	  boost::optional<std::vector<LabelSet>>());
                 //		}
                 //
                 //		// Does the vector have the right size?
@@ -164,16 +164,16 @@ namespace storm {
                 //			newTransitionRewards = newTransRewardsBuilder.build();
                 //		}
                 //
-                //		boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> newChoiceLabels;
+                //		boost::optional<std::vector<LabelSet>> newChoiceLabels;
                 //		if(this->hasChoiceLabeling()) {
                 //
                 //			// Get the choice label sets and move the needed values to the front.
-                //			std::vector<boost::container::flat_set<uint_fast64_t>> newChoice(this->getChoiceLabeling());
+                //			std::vector<LabelSet> newChoice(this->getChoiceLabeling());
                 //			storm::utility::vector::selectVectorValues(newChoice, subSysStates, this->getChoiceLabeling());
                 //
                 //			// Throw away all values after the last state and set the choice label set for s_b as empty.
                 //			newChoice.resize(newStateCount);
-                //			newChoice[newStateCount - 1] = boost::container::flat_set<uint_fast64_t>();
+                //			newChoice[newStateCount - 1] = LabelSet();
                 //
                 //			newChoiceLabels = newChoice;
                 //		}
@@ -182,6 +182,46 @@ namespace storm {
                 //		return storm::models::Dtmc<ValueType>(newMatBuilder.build(), newLabeling, newStateRewards, std::move(newTransitionRewards), newChoiceLabels);
             }
             
+            template<typename ValueType>
+            Dtmc<ValueType>::ConstraintCollector::ConstraintCollector(Dtmc<ValueType> const& dtmc) {
+                process(dtmc);
+            }
+            
+            template<typename ValueType>
+            std::unordered_set<carl::Constraint<ValueType>> const& Dtmc<ValueType>::ConstraintCollector::getWellformedConstraints() const {
+                return this->wellformedConstraintSet;
+            }
+            
+            template<typename ValueType>
+            std::unordered_set<carl::Constraint<ValueType>> const& Dtmc<ValueType>::ConstraintCollector::getGraphPreservingConstraints() const {
+                return this->graphPreservingConstraintSet;
+            }
+            
+            template<typename ValueType>
+            void Dtmc<ValueType>::ConstraintCollector::process(storm::models::sparse::Dtmc<ValueType> const& dtmc) {
+                for(uint_fast64_t state = 0; state < dtmc.getNumberOfStates(); ++state) {
+                    ValueType sum = storm::utility::zero<ValueType>();
+                    for (auto const& transition : dtmc.getRows(state)) {
+                        sum += transition.getValue();
+                        if (!comparator.isConstant(transition.getValue())) {
+                            wellformedConstraintSet.emplace(transition.getValue() - 1, storm::CompareRelation::LEQ);
+                            wellformedConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GEQ);
+                            graphPreservingConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GT);
+                        }
+                    }
+                    STORM_LOG_ASSERT(!comparator.isConstant(sum) || comparator.isOne(sum), "If the sum is a constant, it must be equal to 1.");
+                    if(!comparator.isConstant(sum)) {
+                        wellformedConstraintSet.emplace(sum - 1, storm::CompareRelation::EQ);
+                    }
+                    
+                }
+            }
+            
+            template<typename ValueType>
+            void Dtmc<ValueType>::ConstraintCollector::operator()(storm::models::sparse::Dtmc<ValueType> const& dtmc) {
+                process(dtmc);
+            }
+            
             template <typename ValueType>
             bool Dtmc<ValueType>::checkValidityOfProbabilityMatrix() const {
                 if (this->getTransitionMatrix().getRowCount() != this->getTransitionMatrix().getColumnCount()) {

src/models/sparse/Dtmc.h

@@ -27,7 +27,7 @@ namespace storm {
                      storm::models::sparse::StateLabeling const& stateLabeling,
                      boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                      boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                     boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                     boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -41,7 +41,7 @@ namespace storm {
                 Dtmc(storm::storage::SparseMatrix<ValueType>&& probabilityMatrix, storm::models::sparse::StateLabeling&& stateLabeling,
                      boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                      boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                     boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                     boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 Dtmc(Dtmc<ValueType> const& dtmc) = default;
                 Dtmc& operator=(Dtmc<ValueType> const& dtmc) = default;
@@ -59,6 +59,57 @@ namespace storm {
                  */
                 Dtmc<ValueType> getSubDtmc(storm::storage::BitVector const& states) const;
                 
+                class ConstraintCollector {
+                private:
+                    // A set of constraints that says that the DTMC actually has valid probability distributions in all states.
+                    std::unordered_set<carl::Constraint<ValueType>> wellformedConstraintSet;
+                    
+                    // A set of constraints that makes sure that the underlying graph of the model does not change depending
+                    // on the parameter values.
+                    std::unordered_set<carl::Constraint<ValueType>> graphPreservingConstraintSet;
+                    
+                    // A comparator that is used for
+                    storm::utility::ConstantsComparator<ValueType> comparator;
+                    
+                public:
+                    /*!
+                     * Constructs the a constraint collector for the given DTMC. The constraints are built and ready for
+                     * retrieval after the construction.
+                     *
+                     * @param dtmc The DTMC for which to create the constraints.
+                     */
+                    ConstraintCollector(storm::models::sparse::Dtmc<ValueType> const& dtmc);
+                    
+                    /*!
+                     * Returns the set of wellformed-ness constraints.
+                     *
+                     * @return The set of wellformed-ness constraints.
+                     */
+                    std::unordered_set<carl::Constraint<ValueType>> const&  getWellformedConstraints() const;
+                    
+                    /*!
+                     * Returns the set of graph-preserving constraints.
+                     *
+                     * @return The set of graph-preserving constraints.
+                     */
+                    std::unordered_set<carl::Constraint<ValueType>> const&  getGraphPreservingConstraints() const;
+                    
+                    /*!
+                     * Constructs the constraints for the given DTMC.
+                     *
+                     * @param dtmc The DTMC for which to create the constraints.
+                     */
+                    void process(storm::models::sparse::Dtmc<ValueType> const& dtmc);
+                    
+                    /*!
+                     * Constructs the constraints for the given DTMC by calling the process method.
+                     *
+                     * @param dtmc The DTMC for which to create the constraints.
+                     */
+                    void operator()(storm::models::sparse::Dtmc<ValueType> const& dtmc);
+                    
+                };
+                
             private:
                 /*!
                  * Checks the probability matrix for validity.

src/models/sparse/MarkovAutomaton.cpp

@@ -13,7 +13,7 @@ namespace storm {
                                                         std::vector<ValueType> const& exitRates,
                                                         boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                                         boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                                        boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                                        boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : NondeterministicModel<ValueType>(storm::models::ModelType::MarkovAutomaton, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling), markovianStates(markovianStates), exitRates(exitRates), closed(false) {
                 this->turnRatesToProbabilities();
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "There are transition rewards for nonexistent transitions.");
@@ -26,7 +26,7 @@ namespace storm {
                                                         std::vector<ValueType> const& exitRates,
                                                         boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                                         boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                                        boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                                        boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : NondeterministicModel<ValueType>(storm::models::ModelType::MarkovAutomaton, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)), markovianStates(markovianStates), exitRates(std::move(exitRates)), closed(false) {
                 this->turnRatesToProbabilities();
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "There are transition rewards for nonexistent transitions.");

src/models/sparse/MarkovAutomaton.h

@@ -31,7 +31,7 @@ namespace storm {
                                 std::vector<ValueType> const& exitRates,
                                 boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                 boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -50,7 +50,7 @@ namespace storm {
                                 std::vector<ValueType> const& exitRates,
                                 boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                 boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 MarkovAutomaton(MarkovAutomaton const& other) = default;
                 MarkovAutomaton& operator=(MarkovAutomaton const& other) = default;

src/models/sparse/Mdp.cpp

@@ -11,7 +11,7 @@ namespace storm {
                                 storm::models::sparse::StateLabeling const& stateLabeling,
                                 boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                 boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : NondeterministicModel<ValueType>(storm::models::ModelType::Mdp, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
                 STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid.");
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
@@ -23,20 +23,20 @@ namespace storm {
                                 storm::models::sparse::StateLabeling&& stateLabeling,
                                 boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                 boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : NondeterministicModel<ValueType>(storm::models::ModelType::Mdp, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) {
                 STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid.");
                 STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
             }
             
             template <typename ValueType>
-            Mdp<ValueType> Mdp<ValueType>::restrictChoiceLabels(boost::container::flat_set<uint_fast64_t> const& enabledChoiceLabels) const {
+            Mdp<ValueType> Mdp<ValueType>::restrictChoiceLabels(LabelSet const& enabledChoiceLabels) const {
                 STORM_LOG_THROW(this->hasChoiceLabeling(), storm::exceptions::InvalidArgumentException, "Restriction to label set is impossible for unlabeled model.");
                 
-                std::vector<boost::container::flat_set<uint_fast64_t>> const& choiceLabeling = this->getChoiceLabeling();
+                std::vector<LabelSet> const& choiceLabeling = this->getChoiceLabeling();
                 
                 storm::storage::SparseMatrixBuilder<ValueType> transitionMatrixBuilder(0, this->getTransitionMatrix().getColumnCount(), 0, true, true);
-                std::vector<boost::container::flat_set<uint_fast64_t>> newChoiceLabeling;
+                std::vector<LabelSet> newChoiceLabeling;
                 
                 // Check for each choice of each state, whether the choice labels are fully contained in the given label set.
                 uint_fast64_t currentRow = 0;
@@ -71,7 +71,7 @@ namespace storm {
                 Mdp<ValueType> restrictedMdp(transitionMatrixBuilder.build(), storm::models::sparse::StateLabeling(this->getStateLabeling()),
                                              this->hasStateRewards() ? boost::optional<std::vector<ValueType>>(this->getStateRewardVector()) : boost::optional<std::vector<ValueType>>(),
                                              this->hasTransitionRewards() ? boost::optional<storm::storage::SparseMatrix<ValueType>>(this->getTransitionRewardMatrix()) : boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                             boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>(newChoiceLabeling));
+                                             boost::optional<std::vector<LabelSet>>(newChoiceLabeling));
                 return restrictedMdp;
             }
             

src/models/sparse/Mdp.h

@@ -27,7 +27,7 @@ namespace storm {
                     storm::models::sparse::StateLabeling const& stateLabeling,
                     boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                     boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                    boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -42,7 +42,7 @@ namespace storm {
                     storm::models::sparse::StateLabeling&& stateLabeling,
                     boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                     boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                    boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 Mdp(Mdp const& other) = default;
                 Mdp& operator=(Mdp const& other) = default;
@@ -61,7 +61,7 @@ namespace storm {
                  * and which ones need to be ignored.
                  * @return A restricted version of the current MDP that only uses choice labels from the given set.
                  */
-                Mdp<ValueType> restrictChoiceLabels(boost::container::flat_set<uint_fast64_t> const& enabledChoiceLabels) const;
+                Mdp<ValueType> restrictChoiceLabels(LabelSet const& enabledChoiceLabels) const;
                 
             private:
                 /*!

src/models/sparse/Model.cpp

@@ -12,7 +12,7 @@ namespace storm {
                                     storm::models::sparse::StateLabeling const& stateLabeling,
                                     boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                     boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                    boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : ModelBase(modelType),
             transitionMatrix(transitionMatrix),
             stateLabeling(stateLabeling),
@@ -28,7 +28,7 @@ namespace storm {
                                     storm::models::sparse::StateLabeling&& stateLabeling,
                                     boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                     boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                    boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : ModelBase(modelType),
             transitionMatrix(std::move(transitionMatrix)),
             stateLabeling(std::move(stateLabeling)),
@@ -109,7 +109,7 @@ namespace storm {
             }
             
             template <typename ValueType>
-            std::vector<boost::container::flat_set<uint_fast64_t>> const& Model<ValueType>::getChoiceLabeling() const {
+            std::vector<LabelSet> const& Model<ValueType>::getChoiceLabeling() const {
                 return choiceLabeling.get();
             }
             
@@ -159,7 +159,7 @@ namespace storm {
                     result += getTransitionRewardMatrix().getSizeInBytes();
                 }
                 if (hasChoiceLabeling()) {
-                    result += getChoiceLabeling().size() * sizeof(boost::container::flat_set<uint_fast64_t>);
+                    result += getChoiceLabeling().size() * sizeof(LabelSet);
                 }
                 return result;
             }

src/models/sparse/Model.h

@@ -16,6 +16,9 @@ namespace storm {
     namespace models {
         namespace sparse {
             
+            // The type used for storing a set of labels.
+            typedef boost::container::flat_set<uint_fast64_t> LabelSet;
+            
             /*!
              * Base class for all sparse models.
              */
@@ -45,7 +48,7 @@ namespace storm {
                       storm::models::sparse::StateLabeling const& stateLabeling,
                       boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                       boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                      boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                      boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -62,7 +65,7 @@ namespace storm {
                       storm::models::sparse::StateLabeling&& stateLabeling,
                       boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                       boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                      boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                      boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                                 
                 /*!
                  * Retrieves the backward transition relation of the model, i.e. a set of transitions between states
@@ -175,7 +178,7 @@ namespace storm {
                  *
                  * @return The labels for the choices of the model.
                  */
-                std::vector<boost::container::flat_set<uint_fast64_t>> const& getChoiceLabeling() const;
+                std::vector<LabelSet> const& getChoiceLabeling() const;
                 
                 /*!
                  * Returns the state labeling associated with this model.
@@ -287,7 +290,7 @@ namespace storm {
                 boost::optional<storm::storage::SparseMatrix<ValueType>> transitionRewardMatrix;
                 
                 // If set, a vector representing the labels of choices.
-                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> choiceLabeling;
+                boost::optional<std::vector<LabelSet>> choiceLabeling;
             };
             
         } // namespace sparse

src/models/sparse/NondeterministicModel.cpp

@@ -12,7 +12,7 @@ namespace storm {
                                                                     storm::models::sparse::StateLabeling const& stateLabeling,
                                                                     boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
                                                                     boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix,
-                                                                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+                                                                    boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
             : Model<ValueType>(modelType, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
                 // Intentionally left empty.
             }
@@ -23,7 +23,7 @@ namespace storm {
                                                                     storm::models::sparse::StateLabeling&& stateLabeling,
                                                                     boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
                                                                     boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix,
-                                                                    boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                                                                    boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
             : Model<ValueType>(modelType, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix),
                                std::move(optionalChoiceLabeling)) {
                 // Intentionally left empty.

src/models/sparse/NondeterministicModel.h

@@ -29,7 +29,7 @@ namespace storm {
                                       storm::models::sparse::StateLabeling const& stateLabeling,
                                       boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                       boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                      boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                      boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 /*!
                  * Constructs a model by moving the given data.
@@ -46,7 +46,7 @@ namespace storm {
                                       storm::models::sparse::StateLabeling&& stateLabeling,
                                       boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
                                       boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(),
-                                      boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
+                                      boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
                 
                 NondeterministicModel(NondeterministicModel const& other) = default;
                 NondeterministicModel& operator=(NondeterministicModel const& other) = default;

src/models/sparse/StochasticTwoPlayerGame.cpp

@@ -0,0 +1,41 @@
+#include "src/models/sparse/StochasticTwoPlayerGame.h"
+
+#include "src/adapters/CarlAdapter.h"
+
+namespace storm {
+    namespace models {
+        namespace sparse {
+            
+            template <typename ValueType>
+            StochasticTwoPlayerGame<ValueType>::StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix,
+                                                                        storm::storage::SparseMatrix<ValueType> const& player2Matrix,
+                                                                        storm::models::sparse::StateLabeling const& stateLabeling,
+                                                                        boost::optional<std::vector<ValueType>> const& optionalStateRewardVector,
+                                                                        boost::optional<std::vector<LabelSet>> const& optionalPlayer1ChoiceLabeling,
+                                                                        boost::optional<std::vector<LabelSet>> const& optionalPlayer2ChoiceLabeling)
+            : NondeterministicModel<ValueType>(storm::models::ModelType::S2pg, player2Matrix, stateLabeling, optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>(), optionalPlayer2ChoiceLabeling), player1Matrix(player1Matrix), player1Labels(optionalPlayer1ChoiceLabeling) {
+                // Intentionally left empty.
+            }
+            
+            
+            template <typename ValueType>
+            StochasticTwoPlayerGame<ValueType>::StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type>&& player1Matrix,
+                                                                        storm::storage::SparseMatrix<ValueType>&& player2Matrix,
+                                                                        storm::models::sparse::StateLabeling&& stateLabeling,
+                                                                        boost::optional<std::vector<ValueType>>&& optionalStateRewardVector,
+                                                                        boost::optional<std::vector<LabelSet>>&& optionalPlayer1ChoiceLabeling,
+                                                                        boost::optional<std::vector<LabelSet>>&& optionalPlayer2ChoiceLabeling)
+            : NondeterministicModel<ValueType>(storm::models::ModelType::S2pg, std::move(player2Matrix), std::move(stateLabeling), std::move(optionalStateRewardVector), boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(optionalPlayer2ChoiceLabeling)), player1Matrix(std::move(player1Matrix)), player1Labels(std::move(optionalPlayer1ChoiceLabeling)) {
+                // Intentionally left empty.
+            }
+            
+            template class StochasticTwoPlayerGame<double>;
+            template class StochasticTwoPlayerGame<float>;
+            
+#ifdef STORM_HAVE_CARL
+            template class StochasticTwoPlayerGame<storm::RationalFunction>;
+#endif
+            
+        } // namespace sparse
+    } // namespace models
+} // namespace storm

src/models/sparse/StochasticTwoPlayerGame.h

@@ -0,0 +1,77 @@
+#ifndef STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_
+#define STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_
+
+#include "src/models/sparse/NondeterministicModel.h"
+#include "src/utility/OsDetection.h"
+
+namespace storm {
+    namespace models {
+        namespace sparse {
+            
+            /*!
+             * This class represents a (discrete-time) stochastic two-player game.
+             */
+            template <typename ValueType>
+            class StochasticTwoPlayerGame : public NondeterministicModel<ValueType> {
+            public:
+                /*!
+                 * Constructs a model from the given data.
+                 *
+                 * @param player1Matrix The matrix representing the choices of player 1.
+                 * @param player2Matrix The matrix representing the choices of player 2.
+                 * @param stateLabeling The labeling of the states.
+                 * @param optionalStateRewardVector The reward values associated with the states.
+                 * @param optionalPlayer1ChoiceLabeling A vector that represents the labels associated with the choices of each player 1 state.
+                 * @param optionalPlayer2ChoiceLabeling A vector that represents the labels associated with the choices of each player 2 state.
+                 */
+                StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix,
+                                        storm::storage::SparseMatrix<ValueType> const& player2Matrix,
+                                        storm::models::sparse::StateLabeling const& stateLabeling,
+                                        boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
+                                        boost::optional<std::vector<LabelSet>> const& optionalPlayer1ChoiceLabeling = boost::optional<std::vector<LabelSet>>(),
+                                        boost::optional<std::vector<LabelSet>> const& optionalPlayer2ChoiceLabeling = boost::optional<std::vector<LabelSet>>());
+                
+                /*!
+                 * Constructs a model by moving the given data.
+                 *
+                 * @param player1Matrix The matrix representing the choices of player 1.
+                 * @param player2Matrix The matrix representing the choices of player 2.
+                 * @param stateLabeling The labeling of the states.
+                 * @param optionalStateRewardVector The reward values associated with the states.
+                 * @param optionalPlayer1ChoiceLabeling A vector that represents the labels associated with the choices of each player 1 state.
+                 * @param optionalPlayer2ChoiceLabeling A vector that represents the labels associated with the choices of each player 2 state.
+                 */
+                StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type>&& player1Matrix,
+                                        storm::storage::SparseMatrix<ValueType>&& player2Matrix,
+                                        storm::models::sparse::StateLabeling&& stateLabeling,
+                                        boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(),
+                                        boost::optional<std::vector<LabelSet>>&& optionalPlayer1ChoiceLabeling = boost::optional<std::vector<LabelSet>>(),
+                                        boost::optional<std::vector<LabelSet>>&& optionalPlayer2ChoiceLabeling = boost::optional<std::vector<LabelSet>>());
+                
+                StochasticTwoPlayerGame(StochasticTwoPlayerGame const& other) = default;
+                StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame const& other) = default;
+                
+#ifndef WINDOWS
+                StochasticTwoPlayerGame(StochasticTwoPlayerGame&& other) = default;
+                StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame&& other) = default;
+#endif
+                
+            private:
+                // A matrix that stores the player 1 choices. This matrix contains a row group for each player 1 node. Every
+                // row group contains a row for each choice in that player 1 node. Each such row contains exactly one
+                // (non-zero) entry at a column that indicates the player 2 node this choice leads to (which is essentially
+                // the index of a row group in the matrix for player 2).
+                storm::storage::SparseMatrix<storm::storage::sparse::state_type> player1Matrix;
+                
+                // An (optional) vector of labels attached to the choices of player 1. Each row of the matrix can be equipped
+                // with a set of labels to tag certain choices.
+                boost::optional<std::vector<LabelSet>> player1Labels;
+                
+                // The matrix and labels for player 2 are stored in the superclass.
+            };
+            
+        } // namespace sparse
+    } // namespace models
+} // namespace storm
+
+#endif /* STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_ */

src/models/symbolic/StochasticTwoPlayerGame.cpp

@@ -0,0 +1,42 @@
+#include "src/models/symbolic/StochasticTwoPlayerGame.h"
+
+namespace storm {
+    namespace models {
+        namespace symbolic {
+            
+            template<storm::dd::DdType Type>
+            StochasticTwoPlayerGame<Type>::StochasticTwoPlayerGame(std::shared_ptr<storm::dd::DdManager<Type>> manager,
+                                                                   storm::dd::Bdd<Type> reachableStates,
+                                                                   storm::dd::Bdd<Type> initialStates,
+                                                                   storm::dd::Add<Type> transitionMatrix,
+                                                                   std::set<storm::expressions::Variable> const& rowVariables,
+                                                                   std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> rowExpressionAdapter,
+                                                                   std::set<storm::expressions::Variable> const& columnVariables,
+                                                                   std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> columnExpressionAdapter,
+                                                                   std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs,
+                                                                   std::set<storm::expressions::Variable> const& player1Variables,
+                                                                   std::set<storm::expressions::Variable> const& player2Variables,
+                                                                   std::set<storm::expressions::Variable> const& nondeterminismVariables,
+                                                                   std::map<std::string, storm::expressions::Expression> labelToExpressionMap,
+                                                                   boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector,
+                                                                   boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix)
+                : NondeterministicModel<Type>(storm::models::ModelType::S2pg, manager, reachableStates, initialStates, transitionMatrix, rowVariables, rowExpressionAdapter, columnVariables, columnExpressionAdapter, rowColumnMetaVariablePairs, nondeterminismVariables, labelToExpressionMap, optionalStateRewardVector, optionalTransitionRewardMatrix), player1Variables(player1Variables), player2Variables(player2Variables) {
+                // Intentionally left empty.
+            }
+            
+            template<storm::dd::DdType Type>
+            std::set<storm::expressions::Variable> const& StochasticTwoPlayerGame<Type>::getPlayer1Variables() const {
+                return player1Variables;
+            }
+            
+            template<storm::dd::DdType Type>
+            std::set<storm::expressions::Variable> const& StochasticTwoPlayerGame<Type>::getPlayer2Variables() const {
+                return player2Variables;
+            }
+            
+            // Explicitly instantiate the template class.
+            template class StochasticTwoPlayerGame<storm::dd::DdType::CUDD>;
+            
+        } // namespace symbolic
+    } // namespace models
+} // namespace storm

src/models/symbolic/StochasticTwoPlayerGame.h

@@ -0,0 +1,88 @@
+#ifndef STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_
+#define STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_
+
+#include "src/models/symbolic/NondeterministicModel.h"
+#include "src/utility/OsDetection.h"
+
+namespace storm {
+    namespace models {
+        namespace symbolic {
+            
+            /*!
+             * This class represents a discrete-time stochastic two-player game.
+             */
+            template<storm::dd::DdType Type>
+            class StochasticTwoPlayerGame : public NondeterministicModel<Type> {
+            public:
+                StochasticTwoPlayerGame(StochasticTwoPlayerGame<Type> const& other) = default;
+                StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame<Type> const& other) = default;
+                
+#ifndef WINDOWS
+                StochasticTwoPlayerGame(StochasticTwoPlayerGame<Type>&& other) = default;
+                StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame<Type>&& other) = default;
+#endif
+                
+                /*!
+                 * Constructs a model from the given data.
+                 *
+                 * @param manager The manager responsible for the decision diagrams.
+                 * @param reachableStates A DD representing the reachable states.
+                 * @param initialStates A DD representing the initial states of the model.
+                 * @param transitionMatrix The matrix representing the transitions in the model.
+                 * @param rowVariables The set of row meta variables used in the DDs.
+                 * @param rowExpressionAdapter An object that can be used to translate expressions in terms of the row
+                 * meta variables.
+                 * @param columVariables The set of column meta variables used in the DDs.
+                 * @param columnExpressionAdapter An object that can be used to translate expressions in terms of the
+                 * column meta variables.
+                 * @param rowColumnMetaVariablePairs All pairs of row/column meta variables.
+                 * @param player1Variables The meta variables used to encode the nondeterministic choices of player 1.
+                 * @param player2Variables The meta variables used to encode the nondeterministic choices of player 2.
+                 * @param allNondeterminismVariables The meta variables used to encode the nondeterminism in the model.
+                 * @param labelToExpressionMap A mapping from label names to their defining expressions.
+                 * @param optionalStateRewardVector The reward values associated with the states.
+                 * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model.
+                 */
+                StochasticTwoPlayerGame(std::shared_ptr<storm::dd::DdManager<Type>> manager,
+                                        storm::dd::Bdd<Type> reachableStates,
+                                        storm::dd::Bdd<Type> initialStates,
+                                        storm::dd::Add<Type> transitionMatrix,
+                                        std::set<storm::expressions::Variable> const& rowVariables,
+                                        std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> rowExpressionAdapter,
+                                        std::set<storm::expressions::Variable> const& columnVariables,
+                                        std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> columnExpressionAdapter,
+                                        std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs,
+                                        std::set<storm::expressions::Variable> const& player1Variables,
+                                        std::set<storm::expressions::Variable> const& player2Variables,
+                                        std::set<storm::expressions::Variable> const& allNondeterminismVariables,
+                                        std::map<std::string, storm::expressions::Expression> labelToExpressionMap = std::map<std::string, storm::expressions::Expression>(),
+                                        boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector = boost::optional<storm::dd::Dd<Type>>(),
+                                        boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix = boost::optional<storm::dd::Dd<Type>>());
+
+                /*!
+                 * Retrieeves the set of meta variables used to encode the nondeterministic choices of player 1.
+                 *
+                 * @return The set of meta variables used to encode the nondeterministic choices of player 1.
+                 */
+                std::set<storm::expressions::Variable> const& getPlayer1Variables() const;
+                
+                /*!
+                 * Retrieeves the set of meta variables used to encode the nondeterministic choices of player 2.
+                 *
+                 * @return The set of meta variables used to encode the nondeterministic choices of player 2.
+                 */
+                std::set<storm::expressions::Variable> const& getPlayer2Variables() const;
+                
+            private:
+                // The meta variables used to encode the nondeterministic choices of player 1.
+                std::set<storm::expressions::Variable> player1Variables;
+                
+                // The meta variables used to encode the nondeterministic choices of player 2.
+                std::set<storm::expressions::Variable> player2Variables;
+            };
+            
+        } // namespace symbolic
+    } // namespace models
+} // namespace storm
+
+#endif /* STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_ */

src/parser/ExpressionParser.cpp

@@ -23,7 +23,7 @@ namespace storm {
             }
             minMaxExpression.name("min/max expression");
             
-            identifierExpression = identifier[qi::_val = phoenix::bind(&ExpressionParser::getIdentifierExpression, phoenix::ref(*this), qi::_1)];
+            identifierExpression = identifier[qi::_val = phoenix::bind(&ExpressionParser::getIdentifierExpression, phoenix::ref(*this), qi::_1, allowBacktracking, phoenix::ref(qi::_pass))];
             identifierExpression.name("identifier expression");
             
             literalExpression = trueFalse_[qi::_val = qi::_1] | strict_double[qi::_val = phoenix::bind(&ExpressionParser::createDoubleLiteralExpression, phoenix::ref(*this), qi::_1, qi::_pass)] | qi::int_[qi::_val = phoenix::bind(&ExpressionParser::createIntegerLiteralExpression, phoenix::ref(*this), qi::_1)];
@@ -53,7 +53,7 @@ namespace storm {
             plusExpression.name("plus expression");
             
             if (allowBacktracking) {
-                relativeExpression = plusExpression[qi::_val = qi::_1] > -(relationalOperator_ >> plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
+                relativeExpression = plusExpression[qi::_val = qi::_1] >> -(relationalOperator_ >> plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
             } else {
                 relativeExpression = plusExpression[qi::_val = qi::_1] > -(relationalOperator_ > plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
             }
@@ -70,7 +70,7 @@ namespace storm {
             andExpression.name("and expression");
             
             if (allowBacktracking) {
-                orExpression = andExpression[qi::_val = qi::_1] > *(orOperator_ >> andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
+                orExpression = andExpression[qi::_val = qi::_1] >> *(orOperator_ >> andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
             } else {
                 orExpression = andExpression[qi::_val = qi::_1] > *(orOperator_ > andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
             }
@@ -82,6 +82,24 @@ namespace storm {
             expression %= iteExpression;
             expression.name("expression");
             
+            /*!
+             * Enable this for debugging purposes.
+            debug(expression);
+            debug(iteExpression);
+            debug(orExpression);
+            debug(andExpression);
+            debug(equalityExpression);
+            debug(relativeExpression);
+            debug(plusExpression);
+            debug(multiplicationExpression);
+            debug(powerExpression);
+            debug(unaryExpression);
+            debug(atomicExpression);
+            debug(literalExpression);
+            debug(identifierExpression);
+             */
+            
+            
             // Enable error reporting.
             qi::on_error<qi::fail>(expression, handler(qi::_1, qi::_2, qi::_3, qi::_4));
             qi::on_error<qi::fail>(iteExpression, handler(qi::_1, qi::_2, qi::_3, qi::_4));
@@ -304,11 +322,18 @@ namespace storm {
             return manager.boolean(false);
         }
         
-        storm::expressions::Expression ExpressionParser::getIdentifierExpression(std::string const& identifier) const {
+        storm::expressions::Expression ExpressionParser::getIdentifierExpression(std::string const& identifier, bool allowBacktracking, bool& pass) const {
             if (this->createExpressions) {
                 STORM_LOG_THROW(this->identifiers_ != nullptr, storm::exceptions::WrongFormatException, "Unable to substitute identifier expressions without given mapping.");
                 storm::expressions::Expression const* expression = this->identifiers_->find(identifier);
-                STORM_LOG_THROW(expression != nullptr, storm::exceptions::WrongFormatException, "Parsing error in line " << get_line(qi::_3) << ": Undeclared identifier '" << identifier << "'.");
+                if (expression == nullptr) {
+                    if (allowBacktracking) {
+                        pass = false;
+                        return manager.boolean(false);
+                    } else {
+                        STORM_LOG_THROW(expression != nullptr, storm::exceptions::WrongFormatException, "Parsing error in line " << get_line(qi::_3) << ": Undeclared identifier '" << identifier << "'.");
+                    }
+                }
                 return *expression;
             } else {
                 return manager.boolean(false);

src/parser/ExpressionParser.h

@@ -223,7 +223,7 @@ namespace storm {
             storm::expressions::Expression createIntegerLiteralExpression(int value) const;
             storm::expressions::Expression createMinimumMaximumExpression(storm::expressions::Expression const& e1, storm::expressions::OperatorType const& operatorType, storm::expressions::Expression const& e2) const;
             storm::expressions::Expression createFloorCeilExpression(storm::expressions::OperatorType const& operatorType, storm::expressions::Expression const& e1) const;
-            storm::expressions::Expression getIdentifierExpression(std::string const& identifier) const;
+            storm::expressions::Expression getIdentifierExpression(std::string const& identifier, bool allowBacktracking, bool& pass) const;
             
             bool isValidIdentifier(std::string const& identifier);
             

src/parser/FormulaParser.cpp

@@ -38,7 +38,10 @@ namespace storm {
             booleanLiteralFormula = (qi::lit("true")[qi::_a = true] | qi::lit("false")[qi::_a = false])[qi::_val = phoenix::bind(&FormulaParser::createBooleanLiteralFormula, phoenix::ref(*this), qi::_a)];
             booleanLiteralFormula.name("boolean literal formula");
             
-            atomicStateFormula = booleanLiteralFormula | labelFormula | expressionFormula | (qi::lit("(") > stateFormula > qi::lit(")"));
+            operatorFormula = probabilityOperator | rewardOperator | steadyStateOperator;
+            operatorFormula.name("operator formulas");
+            
+            atomicStateFormula = booleanLiteralFormula | labelFormula | expressionFormula | (qi::lit("(") > stateFormula > qi::lit(")")) | operatorFormula;
             atomicStateFormula.name("atomic state formula");
 
             notStateFormula = (-unaryBooleanOperator_ >> atomicStateFormula)[qi::_val = phoenix::bind(&FormulaParser::createUnaryBooleanStateFormula, phoenix::ref(*this), qi::_2, qi::_1)];
@@ -89,7 +92,7 @@ namespace storm {
             orStateFormula = andStateFormula[qi::_val = qi::_1] >> *(qi::lit("|") >> andStateFormula)[qi::_val = phoenix::bind(&FormulaParser::createBinaryBooleanStateFormula, phoenix::ref(*this), qi::_val, qi::_1, storm::logic::BinaryBooleanStateFormula::OperatorType::Or)];
             orStateFormula.name("or state formula");
             
-            stateFormula = (probabilityOperator | rewardOperator | steadyStateOperator | orStateFormula);
+            stateFormula = (orStateFormula);
             stateFormula.name("state formula");
             
             start = qi::eps > stateFormula >> qi::skip(boost::spirit::ascii::space | qi::lit("//") >> *(qi::char_ - (qi::eol | qi::eoi)))[qi::eps] >> qi::eoi;

src/parser/FormulaParser.h

@@ -132,6 +132,7 @@ namespace storm {
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> pathFormulaWithoutUntil;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> simplePathFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> atomicStateFormula;
+            qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> operatorFormula;
             qi::rule<Iterator, std::string(), Skipper> label;
             
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> andStateFormula;

src/settings/modules/GmmxxEquationSolverSettings.cpp

@@ -7,7 +7,7 @@ namespace storm {
         namespace modules {
             
             const std::string GmmxxEquationSolverSettings::moduleName = "gmm++";
-            const std::string GmmxxEquationSolverSettings::techniqueOptionName = "tech";
+            const std::string GmmxxEquationSolverSettings::techniqueOptionName = "method";
             const std::string GmmxxEquationSolverSettings::preconditionOptionName = "precond";
             const std::string GmmxxEquationSolverSettings::restartOptionName = "restart";
             const std::string GmmxxEquationSolverSettings::maximalIterationsOptionName = "maxiter";
@@ -32,38 +32,38 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build());
             }
             
-            bool GmmxxEquationSolverSettings::isLinearEquationSystemTechniqueSet() const {
+            bool GmmxxEquationSolverSettings::isLinearEquationSystemMethodSet() const {
                 return this->getOption(techniqueOptionName).getHasOptionBeenSet();
             }
             
-            GmmxxEquationSolverSettings::LinearEquationTechnique GmmxxEquationSolverSettings::getLinearEquationSystemTechnique() const {
+            GmmxxEquationSolverSettings::LinearEquationMethod GmmxxEquationSolverSettings::getLinearEquationSystemMethod() const {
                 std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString();
                 if (linearEquationSystemTechniqueAsString == "bicgstab") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab;
                 } else if (linearEquationSystemTechniqueAsString == "qmr") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Qmr;
                 } else if (linearEquationSystemTechniqueAsString == "gmres") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Gmres;
                 } else if (linearEquationSystemTechniqueAsString == "jacobi") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }
             
-            bool GmmxxEquationSolverSettings::isPreconditioningTechniqueSet() const {
+            bool GmmxxEquationSolverSettings::isPreconditioningMethodSet() const {
                 return this->getOption(preconditionOptionName).getHasOptionBeenSet();
             }
             
-            GmmxxEquationSolverSettings::PreconditioningTechnique GmmxxEquationSolverSettings::getPreconditioningTechnique() const {
-                std::string preconditioningTechniqueAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString();
-                if (preconditioningTechniqueAsString == "ilu") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu;
-                } else if (preconditioningTechniqueAsString == "diagonal") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal;
-                } else if (preconditioningTechniqueAsString == "none") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::None;
+            GmmxxEquationSolverSettings::PreconditioningMethod GmmxxEquationSolverSettings::getPreconditioningMethod() const {
+                std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString();
+                if (PreconditioningMethodAsString == "ilu") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::Ilu;
+                } else if (PreconditioningMethodAsString == "diagonal") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal;
+                } else if (PreconditioningMethodAsString == "none") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::None;
                 }
-                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << preconditioningTechniqueAsString << "' selected.");
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected.");
             }
             
             bool GmmxxEquationSolverSettings::isRestartIterationCountSet() const {
@@ -99,7 +99,7 @@ namespace storm {
             }
             
             bool GmmxxEquationSolverSettings::check() const {
-                bool optionsSet = isLinearEquationSystemTechniqueSet() || isPreconditioningTechniqueSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
+                bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
                 
                 STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
                 

src/settings/modules/GmmxxEquationSolverSettings.h

@@ -12,11 +12,11 @@ namespace storm {
              */
             class GmmxxEquationSolverSettings : public ModuleSettings {
             public:
-                // An enumeration of all available techniques for solving linear equations.
-                enum class LinearEquationTechnique { Bicgstab, Qmr, Gmres, Jacobi };
+                // An enumeration of all available methods for solving linear equations.
+                enum class LinearEquationMethod { Bicgstab, Qmr, Gmres, Jacobi };
 
-                // An enumeration of all available preconditioning techniques.
-                enum class PreconditioningTechnique { Ilu, Diagonal, None };
+                // An enumeration of all available preconditioning methods.
+                enum class PreconditioningMethod { Ilu, Diagonal, None };
                 
                 // An enumeration of all available convergence criteria.
                 enum class ConvergenceCriterion { Absolute, Relative };
@@ -29,32 +29,32 @@ namespace storm {
                 GmmxxEquationSolverSettings(storm::settings::SettingsManager& settingsManager);
                 
                 /*!
-                 * Retrieves whether the linear equation system technique has been set.
+                 * Retrieves whether the linear equation system method has been set.
                  *
-                 * @return True iff the linear equation system technique has been set.
+                 * @return True iff the linear equation system method has been set.
                  */
-                bool isLinearEquationSystemTechniqueSet() const;
+                bool isLinearEquationSystemMethodSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for solving systems of linear equations.
+                 * Retrieves the method that is to be used for solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                LinearEquationTechnique getLinearEquationSystemTechnique() const;
+                LinearEquationMethod getLinearEquationSystemMethod() const;
                 
                 /*!
-                 * Retrieves whether the preconditioning technique has been set.
+                 * Retrieves whether the preconditioning method has been set.
                  *
-                 * @return True iff the preconditioning technique has been set.
+                 * @return True iff the preconditioning method has been set.
                  */
-                bool isPreconditioningTechniqueSet() const;
+                bool isPreconditioningMethodSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for preconditioning solving systems of linear equations.
+                 * Retrieves the method that is to be used for preconditioning solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                PreconditioningTechnique getPreconditioningTechnique() const;
+                PreconditioningMethod getPreconditioningMethod() const;
                 
                 /*!
                  * Retrieves whether the restart iteration count has been set.
@@ -64,7 +64,7 @@ namespace storm {
                 bool isRestartIterationCountSet() const;
                 
                 /*!
-                 * Retrieves the number of iterations after which restarted techniques are to be restarted.
+                 * Retrieves the number of iterations after which restarted methods are to be restarted.
                  *
                  * @return The number of iterations after which to restart.
                  */

src/settings/modules/NativeEquationSolverSettings.cpp

@@ -7,20 +7,23 @@ namespace storm {
         namespace modules {
             
             const std::string NativeEquationSolverSettings::moduleName = "native";
-            const std::string NativeEquationSolverSettings::techniqueOptionName = "tech";
+            const std::string NativeEquationSolverSettings::techniqueOptionName = "method";
+            const std::string NativeEquationSolverSettings::omegaOptionName = "soromega";
             const std::string NativeEquationSolverSettings::maximalIterationsOptionName = "maxiter";
             const std::string NativeEquationSolverSettings::maximalIterationsOptionShortName = "i";
             const std::string NativeEquationSolverSettings::precisionOptionName = "precision";
             const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute";
             
             NativeEquationSolverSettings::NativeEquationSolverSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) {
-                std::vector<std::string> methods = { "jacobi" };
+                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor" };
                 this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine. Available are: { jacobi }.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("jacobi").build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
                 
+                this->addOption(storm::settings::OptionBuilder(moduleName, omegaOptionName, false, "The omega used for SOR.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The value of the SOR parameter.").setDefaultValueDouble(0.9).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+                
                 this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build());
             }
             
@@ -28,10 +31,14 @@ namespace storm {
                 return this->getOption(techniqueOptionName).getHasOptionBeenSet();
             }
             
-            NativeEquationSolverSettings::LinearEquationTechnique NativeEquationSolverSettings::getLinearEquationSystemTechnique() const {
+            NativeEquationSolverSettings::LinearEquationMethod NativeEquationSolverSettings::getLinearEquationSystemMethod() const {
                 std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString();
                 if (linearEquationSystemTechniqueAsString == "jacobi") {
-                    return NativeEquationSolverSettings::LinearEquationTechnique::Jacobi;
+                    return NativeEquationSolverSettings::LinearEquationMethod::Jacobi;
+                } else if (linearEquationSystemTechniqueAsString == "gaussseidel") {
+                    return NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel;
+                } else if (linearEquationSystemTechniqueAsString == "sor") {
+                    return NativeEquationSolverSettings::LinearEquationMethod::SOR;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }
@@ -51,6 +58,10 @@ namespace storm {
             double NativeEquationSolverSettings::getPrecision() const {
                 return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
             }
+
+            double NativeEquationSolverSettings::getOmega() const {
+                return this->getOption(omegaOptionName).getArgumentByName("value").getValueAsDouble();
+            }
             
             bool NativeEquationSolverSettings::isConvergenceCriterionSet() const {
                 return this->getOption(absoluteOptionName).getHasOptionBeenSet();

src/settings/modules/NativeEquationSolverSettings.h

@@ -12,8 +12,8 @@ namespace storm {
              */
             class NativeEquationSolverSettings : public ModuleSettings {
             public:
-                // An enumeration of all available techniques for solving linear equations.
-                enum class LinearEquationTechnique { Jacobi };
+                // An enumeration of all available methods for solving linear equations.
+                enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR };
                 
                 // An enumeration of all available convergence criteria.
                 enum class ConvergenceCriterion { Absolute, Relative };
@@ -33,11 +33,11 @@ namespace storm {
                 bool isLinearEquationSystemTechniqueSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for solving systems of linear equations.
+                 * Retrieves the method that is to be used for solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                LinearEquationTechnique getLinearEquationSystemTechnique() const;
+                LinearEquationMethod getLinearEquationSystemMethod() const;
                 
                 /*!
                  * Retrieves whether the maximal iteration count has been set.
@@ -67,6 +67,13 @@ namespace storm {
                  */
                 double getPrecision() const;
                 
+                /*!
+                 * Retrieves the value of omega to be used for the SOR method.
+                 *
+                 * @return The value of omega to be used.
+                 */
+                double getOmega() const;
+                
                 /*!
                  * Retrieves whether the convergence criterion has been set.
                  *
@@ -89,6 +96,7 @@ namespace storm {
             private:
                 // Define the string names of the options as constants.
                 static const std::string techniqueOptionName;
+                static const std::string omegaOptionName;
                 static const std::string maximalIterationsOptionName;
                 static const std::string maximalIterationsOptionShortName;
                 static const std::string precisionOptionName;

src/solver/GmmxxLinearEquationSolver.cpp

@@ -31,31 +31,31 @@ namespace storm {
             restart = settings.getRestartIterationCount();
             
             // Determine the method to be used.
-            storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique();
-            if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab) {
+            storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod();
+            if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab) {
                 method = SolutionMethod::Bicgstab;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Qmr) {
                 method = SolutionMethod::Qmr;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Gmres) {
                 method = SolutionMethod::Gmres;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi) {
                 method = SolutionMethod::Jacobi;
             }
             
             // Check which preconditioner to use.
-            storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique preconditionAsSetting = settings.getPreconditioningTechnique();
-            if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu) {
+            storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod();
+            if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Ilu) {
                 preconditioner = Preconditioner::Ilu;
-            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal) {
+            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal) {
                 preconditioner = Preconditioner::Diagonal;
-            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::None) {
+            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::None) {
                 preconditioner = Preconditioner::None;
             }
         }
         
         template<typename ValueType>
         void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'.");
+            LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations).");
             if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) {
                 LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
             }

src/solver/NativeLinearEquationSolver.cpp

@@ -15,7 +15,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method) : A(A), method(method) {
             // Get the settings object to customize linear solving.
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
             
@@ -23,61 +23,94 @@ namespace storm {
             maximalNumberOfIterations = settings.getMaximalIterationCount();
             precision = settings.getPrecision();
             relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
-            
-            // Determine the method to be used.
-            storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique();
-            if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique::Jacobi) {
-                method = SolutionMethod::Jacobi;
-            }
         }
                 
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            // Get a Jacobi decomposition of the matrix A.
-            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
-            
-            // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
-            bool multiplyResultProvided = true;
-            std::vector<ValueType>* nextX = multiplyResult;
-            if (nextX == nullptr) {
-                nextX = new std::vector<ValueType>(x.size());
-                multiplyResultProvided = false;
-            }
-            std::vector<ValueType> const* copyX = nextX;
-            std::vector<ValueType>* currentX = &x;
-            
-            // Target vector for precision calculation.
-            std::vector<ValueType> tmpX(x.size());
-            
-            // Set up additional environment variables.
-            uint_fast64_t iterationCount = 0;
-            bool converged = false;
-            
-            while (!converged && iterationCount < maximalNumberOfIterations) {
-                // Compute D^-1 * (b - LU * x) and store result in nextX.
-                jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
-                storm::utility::vector::subtractVectors(b, tmpX, tmpX);
-                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
+            if (method == SolutionMethod::SOR || method == SolutionMethod::GaussSeidel) {
+                // Define the omega used for SOR.
+                ValueType omega = method == SolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one<ValueType>();
                 
-                // Swap the two pointers as a preparation for the next iteration.
-                std::swap(nextX, currentX);
+                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+                bool tmpXProvided = true;
+                std::vector<ValueType>* tmpX = multiplyResult;
+                if (multiplyResult == nullptr) {
+                    tmpX = new std::vector<ValueType>(x);
+                    tmpXProvided = false;
+                } else {
+                    *tmpX = x;
+                }
                 
-                // Now check if the process already converged within our precision.
-                converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative);
+                // Set up additional environment variables.
+                uint_fast64_t iterationCount = 0;
+                bool converged = false;
                 
-                // Increase iteration count so we can abort if convergence is too slow.
-                ++iterationCount;
-            }
-            
-            // If the last iteration did not write to the original x we have to swap the contents, because the
-            // output has to be written to the input parameter x.
-            if (currentX == copyX) {
-                std::swap(x, *currentX);
-            }
-            
-            // If the vector for the temporary multiplication result was not provided, we need to delete it.
-            if (!multiplyResultProvided) {
-                delete copyX;
+                while (!converged && iterationCount < maximalNumberOfIterations) {
+                    A.performSuccessiveOverRelaxationStep(omega, x, b);
+                    
+                    // Now check if the process already converged within our precision.
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(precision), relative);
+                    
+                    // If we did not yet converge, we need to copy the contents of x to *tmpX.
+                    if (!converged) {
+                        *tmpX = x;
+                    }
+                    
+                    // Increase iteration count so we can abort if convergence is too slow.
+                    ++iterationCount;
+                }
+                
+                // If the vector for the temporary multiplication result was not provided, we need to delete it.
+                if (!tmpXProvided) {
+                    delete tmpX;
+                }
+            } else {
+                // Get a Jacobi decomposition of the matrix A.
+                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+                
+                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+                bool multiplyResultProvided = true;
+                std::vector<ValueType>* nextX = multiplyResult;
+                if (nextX == nullptr) {
+                    nextX = new std::vector<ValueType>(x.size());
+                    multiplyResultProvided = false;
+                }
+                std::vector<ValueType> const* copyX = nextX;
+                std::vector<ValueType>* currentX = &x;
+                
+                // Target vector for precision calculation.
+                std::vector<ValueType> tmpX(x.size());
+                
+                // Set up additional environment variables.
+                uint_fast64_t iterationCount = 0;
+                bool converged = false;
+                
+                while (!converged && iterationCount < maximalNumberOfIterations) {
+                    // Compute D^-1 * (b - LU * x) and store result in nextX.
+                    jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
+                    storm::utility::vector::subtractVectors(b, tmpX, tmpX);
+                    storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
+                    
+                    // Swap the two pointers as a preparation for the next iteration.
+                    std::swap(nextX, currentX);
+                    
+                    // Now check if the process already converged within our precision.
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative);
+                    
+                    // Increase iteration count so we can abort if convergence is too slow.
+                    ++iterationCount;
+                }
+                
+                // If the last iteration did not write to the original x we have to swap the contents, because the
+                // output has to be written to the input parameter x.
+                if (currentX == copyX) {
+                    std::swap(x, *currentX);
+                }
+                
+                // If the vector for the temporary multiplication result was not provided, we need to delete it.
+                if (!multiplyResultProvided) {
+                    delete copyX;
+                }
             }
         }
         

src/solver/NativeLinearEquationSolver.h

@@ -14,15 +14,16 @@ namespace storm {
         public:
             // An enumeration specifying the available solution methods.
             enum class SolutionMethod {
-                Jacobi
+                Jacobi, GaussSeidel, SOR
             };
             
             /*!
              * Constructs a linear equation solver with parameters being set according to the settings object.
              *
              * @param A The matrix defining the coefficients of the linear equation system.
+             * @param method The method to use for solving linear equations.
              */
-            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method = SolutionMethod::Jacobi);
 
             /*!
              * Constructs a linear equation solver with the given parameters.

src/solver/SymbolicGameSolver.cpp

@@ -0,0 +1,69 @@
+#include "src/solver/SymbolicGameSolver.h"
+
+#include "src/storage/dd/CuddBdd.h"
+#include "src/storage/dd/CuddAdd.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType Type>
+        SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables) {
+            // Get the settings object to customize solving.
+            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
+            storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            precision = settings.getPrecision();
+            relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
+        }
+        
+        template<storm::dd::DdType Type>
+        SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType Type>
+        storm::dd::Add<Type> SymbolicGameSolver<Type>::solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const {
+            // Set up the environment.
+            storm::dd::Add<Type> xCopy = x;
+            uint_fast64_t iterations = 0;
+            bool converged = false;
+            
+            while (!converged && iterations < maximalNumberOfIterations) {
+                // Compute tmp = A * x + b
+                storm::dd::Add<Type> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<Type> tmp = this->gameMatrix.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
+                tmp += b;
+                
+                // Now abstract from player 2 and player 1 variables.
+                if (player2Min) {
+                    tmp = tmp.minAbstract(this->player2Variables);
+                } else {
+                    tmp = tmp.maxAbstract(this->player2Variables);
+                }
+                
+                if (player1Min) {
+                    tmp = tmp.minAbstract(this->player1Variables);
+                } else {
+                    tmp = tmp.maxAbstract(this->player1Variables);
+                }
+                
+                // Now check if the process already converged within our precision.
+                converged = xCopy.equalModuloPrecision(tmp, precision, relative);
+                
+                // If the method did not converge yet, we prepare the x vector for the next iteration.
+                if (!converged) {
+                    xCopy = tmp;
+                }
+                
+                ++iterations;
+            }
+            
+            return xCopy;
+        }
+        
+        template class SymbolicGameSolver<storm::dd::DdType::CUDD>;
+        
+    }
+}

src/solver/SymbolicGameSolver.h

@@ -0,0 +1,97 @@
+#ifndef STORM_SOLVER_SYMBOLICGAMESOLVER_H_
+#define STORM_SOLVER_SYMBOLICGAMESOLVER_H_
+
+#include "src/storage/expressions/Variable.h"
+
+#include "src/storage/dd/Bdd.h"
+#include "src/storage/dd/Add.h"
+
+namespace storm {
+    namespace solver {
+        
+        /*!
+         * An interface that represents an abstract symbolic game solver.
+         */
+        template<storm::dd::DdType Type>
+        class SymbolicGameSolver {
+        public:
+            /*!
+             * Constructs a symbolic game solver with the given meta variable sets and pairs.
+             *
+             * @param gameMatrix The matrix defining the coefficients of the game.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             * @param player1Variables The meta variables used to encode the player 1 choices.
+             * @param player2Variables The meta variables used to encode the player 2 choices.
+             */
+            SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables);
+            
+            /*!
+             * Constructs a symbolic game solver with the given meta variable sets and pairs.
+             *
+             * @param gameMatrix The matrix defining the coefficients of the game.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             * @param player1Variables The meta variables used to encode the player 1 choices.
+             * @param player2Variables The meta variables used to encode the player 2 choices.
+             * @param precision The precision to achieve.
+             * @param maximalNumberOfIterations The maximal number of iterations to perform when solving a linear
+             * equation system iteratively.
+             * @param relative Sets whether or not to use a relativ stopping criterion rather than an absolute one.
+             */
+            SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables, double precision, uint_fast64_t maximalNumberOfIterations, bool relative);
+            
+            /*!
+             * Solves the equation system x = min/max(A*x + b) given by the parameters. Note that the matrix A has
+             * to be given upon construction time of the solver object.
+             *
+             * @param player1Min A flag indicating whether player 1 wants to minimize the result.
+             * @param player2Min A flag indicating whether player 1 wants to minimize the result.
+             * @param x The initial guess of the solution.
+             * @param b The vector to add after matrix-vector multiplication.
+             * @return The solution vector.
+             */
+            virtual storm::dd::Add<Type> solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const;
+        
+        protected:
+            // The matrix defining the coefficients of the linear equation system.
+            storm::dd::Add<Type> const& gameMatrix;
+            
+            // A BDD characterizing all rows of the equation system.
+            storm::dd::Bdd<Type> const& allRows;
+            
+            // The row variables.
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            
+            // The column variables.
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            // The pairs of meta variables used for renaming.
+            std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs;
+            
+            // The player 1 variables.
+            std::set<storm::expressions::Variable> player1Variables;
+            
+            // The player 2 variables.
+            std::set<storm::expressions::Variable> player2Variables;
+            
+            // The precision to achive.
+            double precision;
+            
+            // The maximal number of iterations to perform.
+            uint_fast64_t maximalNumberOfIterations;
+            
+            // A flag indicating whether a relative or an absolute stopping criterion is to be used.
+            bool relative;
+        };
+        
+    } // namespace solver
+} // namespace storm
+
+#endif /* STORM_SOLVER_SYMBOLICGAMESOLVER_H_ */

src/solver/SymbolicLinearEquationSolver.cpp

@@ -0,0 +1,85 @@
+#include "src/solver/SymbolicLinearEquationSolver.h"
+
+#include "src/storage/dd/CuddDdManager.h"
+#include "src/storage/dd/CuddAdd.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
+            // Get the settings object to customize solving.
+            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
+            storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            precision = settings.getPrecision();
+            relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType>  SymbolicLinearEquationSolver<DdType, ValueType>::solveEquationSystem(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const {
+            // Start by computing the Jacobi decomposition of the matrix A.
+            storm::dd::Add<DdType> diagonal = x.getDdManager()->getAddOne();
+            for (auto const& pair : rowColumnMetaVariablePairs) {
+                diagonal *= x.getDdManager()->getIdentity(pair.first).equals(x.getDdManager()->getIdentity(pair.second));
+                diagonal *= x.getDdManager()->getRange(pair.first).toAdd() * x.getDdManager()->getRange(pair.second).toAdd();
+            }
+            diagonal *= allRows.toAdd();
+            
+            storm::dd::Add<DdType> lu = diagonal.ite(this->A.getDdManager()->getAddZero(), this->A);
+            storm::dd::Add<DdType> dinv = diagonal / (diagonal * this->A);
+            
+            // Set up additional environment variables.
+            storm::dd::Add<DdType> xCopy = x;
+            uint_fast64_t iterationCount = 0;
+            bool converged = false;
+            
+            while (!converged && iterationCount < maximalNumberOfIterations) {
+                storm::dd::Add<DdType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<DdType> tmp = lu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
+                tmp = b - tmp;
+                tmp = tmp.swapVariables(this->rowColumnMetaVariablePairs);
+                tmp = dinv.multiplyMatrix(tmp, this->columnMetaVariables);
+                
+                // Now check if the process already converged within our precision.
+                converged = xCopy.equalModuloPrecision(tmp, precision, relative);
+                
+                // If the method did not converge yet, we prepare the x vector for the next iteration.
+                if (!converged) {
+                    xCopy = tmp;
+                }
+                
+                // Increase iteration count so we can abort if convergence is too slow.
+                ++iterationCount;
+            }
+            
+            return xCopy;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> SymbolicLinearEquationSolver<DdType, ValueType>::performMatrixVectorMultiplication(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b, uint_fast64_t n) const {
+            storm::dd::Add<DdType> xCopy = x;
+            
+            // Perform matrix-vector multiplication while the bound is met.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                xCopy = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                xCopy = this->A.multiplyMatrix(xCopy, this->columnMetaVariables);
+                if (b != nullptr) {
+                    xCopy += *b;
+                }
+            }
+            
+            return xCopy;
+        }
+        
+        template class SymbolicLinearEquationSolver<storm::dd::DdType::CUDD, double>;
+        
+    }
+}

src/solver/SymbolicLinearEquationSolver.h

@@ -0,0 +1,104 @@
+#ifndef STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
+
+#include <memory>
+#include <set>
+#include <boost/variant.hpp>
+
+#include "src/storage/expressions/Variable.h"
+#include "src/storage/dd/Bdd.h"
+#include "src/storage/dd/Add.h"
+#include "src/storage/dd/Odd.h"
+
+namespace storm {
+    namespace solver {
+        /*!
+         * An interface that represents an abstract symbolic linear equation solver. In addition to solving a system of
+         * linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
+         */
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicLinearEquationSolver {
+        public:
+            /*!
+             * Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
+             * @param diagonal An ADD characterizing the elements on the diagonal of the matrix.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             */
+            SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+            
+            /*!
+             * Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             * @param precision The precision to achieve.
+             * @param maximalNumberOfIterations The maximal number of iterations to perform when solving a linear
+             * equation system iteratively.
+             * @param relative Sets whether or not to use a relativ stopping criterion rather than an absolute one.
+             */
+            SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative);
+            
+            /*!
+             * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
+             * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
+             * to be given upon construction time of the solver object.
+             *
+             * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
+             * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
+             * @return The solution of the equation system.
+             */
+            virtual storm::dd::Add<DdType> solveEquationSystem(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const;
+            
+            /*!
+             * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
+             * performing the necessary multiplications, the result is written to the input vector x. Note that the
+             * matrix A has to be given upon construction time of the solver object.
+             *
+             * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
+             * to the number of rows of A.
+             * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @return The solution of the equation system.
+             */
+            virtual storm::dd::Add<DdType> performMatrixVectorMultiplication(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b = nullptr, uint_fast64_t n = 1) const;
+            
+        protected:
+            // The matrix defining the coefficients of the linear equation system.
+            storm::dd::Add<DdType> const& A;
+            
+            // A BDD characterizing all rows of the equation system.
+            storm::dd::Bdd<DdType> const& allRows;
+            
+            // The row variables.
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            
+            // The column variables.
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            // The pairs of meta variables used for renaming.
+            std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs;
+            
+            // The precision to achive.
+            double precision;
+            
+            // The maximal number of iterations to perform.
+            uint_fast64_t maximalNumberOfIterations;
+            
+            // A flag indicating whether a relative or an absolute stopping criterion is to be used.
+            bool relative;
+        };
+        
+    } // namespace solver
+} // namespace storm
+
+#endif /* STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_ */

src/solver/Z3SmtSolver.h

@@ -75,7 +75,7 @@ namespace storm {
              * @param model The Z3 model to convert.
              * @return The valuation of variables corresponding to the given model.
              */
-			storm::expressions::SimpleValuation convertZ3ModelToValuation(z3::model const& model);
+            storm::expressions::SimpleValuation convertZ3ModelToValuation(z3::model const& model);
 
             // The context used by the solver.
             std::unique_ptr<z3::context> context;
@@ -87,10 +87,10 @@ namespace storm {
             std::unique_ptr<storm::adapters::Z3ExpressionAdapter> expressionAdapter;
 
             // A flag storing whether the last call to a check method provided aussumptions.
-			bool lastCheckAssumptions;
+            bool lastCheckAssumptions;
             
             // The last result that was returned by any of the check methods.
-			CheckResult lastResult;
+            CheckResult lastResult;
 #endif
 		};
 	}

src/storage/SparseMatrix.cpp

@@ -50,7 +50,7 @@ namespace storm {
         void MatrixEntry<IndexType, ValueType>::setValue(ValueType const& value) {
             this->entry.second = value;
         }
-
+        
         template<typename IndexType, typename ValueType>
         std::pair<IndexType, ValueType> const& MatrixEntry<IndexType, ValueType>::getColumnValuePair() const {
             return this->entry;
@@ -82,13 +82,27 @@ namespace storm {
             rowIndications.push_back(0);
         }
         
+        template<typename ValueType>
+        SparseMatrixBuilder<ValueType>::SparseMatrixBuilder(SparseMatrix<ValueType>&& matrix) :  initialRowCountSet(false), initialRowCount(0), initialColumnCountSet(false), initialColumnCount(0), initialEntryCountSet(false), initialEntryCount(0), forceInitialDimensions(false), hasCustomRowGrouping(matrix.nontrivialRowGrouping), initialRowGroupCountSet(false), initialRowGroupCount(0), rowGroupIndices(), columnsAndValues(std::move(matrix.columnsAndValues)), rowIndications(std::move(matrix.rowIndications)), currentEntryCount(matrix.entryCount), lastRow(matrix.rowCount - 1), lastColumn(columnsAndValues.back().getColumn()), highestColumn(matrix.getColumnCount() - 1), currentRowGroup() {
+            
+            // If the matrix has a custom row grouping, we move it and remove the last element to make it 'open' again.
+            if (hasCustomRowGrouping) {
+                rowGroupIndices = std::move(matrix.rowGroupIndices);
+                rowGroupIndices.pop_back();
+                currentRowGroup = rowGroupIndices.size() - 1;
+            }
+            
+            // Likewise, we need to 'open' the row indications again.
+            rowIndications.pop_back();
+        }
+        
         template<typename ValueType>
         void SparseMatrixBuilder<ValueType>::addNextValue(index_type row, index_type column, ValueType const& value) {
             // Check that we did not move backwards wrt. the row.
             if (row < lastRow) {
                 throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in row " << row << ", but an element in row " << lastRow << " has already been added.";
             }
-
+            
             // Check that we did not move backwards wrt. to column.
             if (row == lastRow && column < lastColumn) {
                 throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in column " << column << " in row " << row << ", but an element in column " << lastColumn << " has already been added in that row.";
@@ -103,14 +117,14 @@ namespace storm {
                 
                 lastRow = row;
             }
-                        
+            
             lastColumn = column;
             
             // Finally, set the element and increase the current size.
             columnsAndValues.emplace_back(column, value);
             highestColumn = std::max(highestColumn, column);
             ++currentEntryCount;
-
+            
             // In case we did not expect this value, we throw an exception.
             if (forceInitialDimensions) {
                 STORM_LOG_THROW(!initialRowCountSet || lastRow < initialRowCount, storm::exceptions::OutOfRangeException, "Cannot insert value at illegal row " << lastRow << ".");
@@ -145,14 +159,14 @@ namespace storm {
             // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for
             // the first and last row.
             rowIndications.push_back(currentEntryCount);
-
+            
             uint_fast64_t columnCount = highestColumn + 1;
             if (initialColumnCountSet && forceInitialDimensions) {
                 STORM_LOG_THROW(columnCount <= initialColumnCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialColumnCount << " columns, but got " << columnCount << ".");
                 columnCount = std::max(columnCount, initialColumnCount);
             }
             columnCount = std::max(columnCount, overriddenColumnCount);
-
+            
             uint_fast64_t entryCount = currentEntryCount;
             if (initialEntryCountSet && forceInitialDimensions) {
                 STORM_LOG_THROW(entryCount == initialEntryCount, storm::exceptions::InvalidStateException, "Expected " << initialEntryCount << " entries, but got " << entryCount << ".");
@@ -170,13 +184,23 @@ namespace storm {
                     rowGroupCount = std::max(rowGroupCount, initialRowGroupCount);
                 }
                 rowGroupCount = std::max(rowGroupCount, overriddenRowGroupCount);
-
+                
                 for (index_type i = currentRowGroup; i <= rowGroupCount; ++i) {
                     rowGroupIndices.push_back(rowCount);
                 }
             }
-
-            return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices));
+            
+            return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), hasCustomRowGrouping);
+        }
+        
+        template<typename ValueType>
+        typename SparseMatrixBuilder<ValueType>::index_type SparseMatrixBuilder<ValueType>::getLastRow() const {
+            return lastRow;
+        }
+        
+        template<typename ValueType>
+        typename SparseMatrixBuilder<ValueType>::index_type SparseMatrixBuilder<ValueType>::getLastColumn() const {
+            return lastColumn;
         }
         
         template<typename ValueType>
@@ -220,17 +244,24 @@ namespace storm {
         }
         
         template<typename ValueType>
-        SparseMatrix<ValueType>::SparseMatrix() : rowCount(0), columnCount(0), entryCount(0), nonzeroEntryCount(0), columnsAndValues(), rowIndications(), rowGroupIndices() {
+        SparseMatrix<ValueType>::SparseMatrix() : rowCount(0), columnCount(0), entryCount(0), nonzeroEntryCount(0), columnsAndValues(), rowIndications(), nontrivialRowGrouping(false), rowGroupIndices() {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType> const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), rowGroupIndices(other.rowGroupIndices) {
+        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType> const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(other.rowGroupIndices) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType>&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), rowGroupIndices(std::move(other.rowGroupIndices)) {
+        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements) {
+            storm::storage::BitVector rowConstraint(other.getRowCount(), true);
+            storm::storage::BitVector columnConstraint(other.getColumnCount(), true);
+            *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements);
+        }
+        
+        template<typename ValueType>
+        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType>&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(std::move(other.rowGroupIndices)) {
             // Now update the source matrix
             other.rowCount = 0;
             other.columnCount = 0;
@@ -238,15 +269,15 @@ namespace storm {
         }
         
         template<typename ValueType>
-        SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) {
-			this->updateNonzeroEntryCount();
+        SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) {
+            this->updateNonzeroEntryCount();
         }
         
         template<typename ValueType>
-        SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) {
-			this->updateNonzeroEntryCount();
+        SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(std::move(rowGroupIndices)) {
+            this->updateNonzeroEntryCount();
         }
-                
+        
         template<typename ValueType>
         SparseMatrix<ValueType>& SparseMatrix<ValueType>::operator=(SparseMatrix<ValueType> const& other) {
             // Only perform assignment if source and target are not the same.
@@ -259,6 +290,7 @@ namespace storm {
                 columnsAndValues = other.columnsAndValues;
                 rowIndications = other.rowIndications;
                 rowGroupIndices = other.rowGroupIndices;
+                nontrivialRowGrouping = other.nontrivialRowGrouping;
             }
             
             return *this;
@@ -276,6 +308,7 @@ namespace storm {
                 columnsAndValues = std::move(other.columnsAndValues);
                 rowIndications = std::move(other.rowIndications);
                 rowGroupIndices = std::move(other.rowGroupIndices);
+                nontrivialRowGrouping = other.nontrivialRowGrouping;
             }
             
             return *this;
@@ -288,7 +321,7 @@ namespace storm {
             }
             
             bool equalityResult = true;
-
+            
             equalityResult &= this->getRowCount() == other.getRowCount();
             if (!equalityResult) {
                 return false;
@@ -346,35 +379,35 @@ namespace storm {
             return entryCount;
         }
         
-		template<typename T>
-		uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const {
-			uint_fast64_t result = 0;
-			for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) {
-				result += (this->rowIndications[row + 1] - this->rowIndications[row]);
-			}
-			return result;
-		}
-
+        template<typename T>
+        uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const {
+            uint_fast64_t result = 0;
+            for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) {
+                result += (this->rowIndications[row + 1] - this->rowIndications[row]);
+            }
+            return result;
+        }
+        
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getNonzeroEntryCount() const {
             return nonzeroEntryCount;
-		}
-
-		template<typename ValueType>
-		void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
-			//SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
-			this->nonzeroEntryCount = 0;
-			for (auto const& element : *this) {
-				if (element.getValue() != storm::utility::zero<ValueType>()) {
-					++this->nonzeroEntryCount;
-				}
-			}
-		}
-
-		template<typename ValueType>
-		void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) {
-			this->nonzeroEntryCount += difference;
-		}
+        }
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
+            //SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
+            this->nonzeroEntryCount = 0;
+            for (auto const& element : *this) {
+                if (element.getValue() != storm::utility::zero<ValueType>()) {
+                    ++this->nonzeroEntryCount;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) {
+            this->nonzeroEntryCount += difference;
+        }
         
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getRowGroupCount() const {
@@ -390,7 +423,7 @@ namespace storm {
         std::vector<typename SparseMatrix<ValueType>::index_type> const& SparseMatrix<ValueType>::getRowGroupIndices() const {
             return rowGroupIndices;
         }
-
+        
         template<typename ValueType>
         void SparseMatrix<ValueType>::makeRowsAbsorbing(storm::storage::BitVector const& rows) {
             for (auto row : rows) {
@@ -424,7 +457,7 @@ namespace storm {
             columnValuePtr->setValue(storm::utility::one<ValueType>());
             ++columnValuePtr;
             for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) {
-				++this->nonzeroEntryCount;
+                ++this->nonzeroEntryCount;
                 columnValuePtr->setColumn(0);
                 columnValuePtr->setValue(storm::utility::zero<ValueType>());
             }
@@ -477,12 +510,47 @@ namespace storm {
                 return getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements);
             }
         }
-                
+        
         template<typename ValueType>
         SparseMatrix<ValueType> SparseMatrix<ValueType>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<index_type> const& rowGroupIndices, bool insertDiagonalEntries) const {
-            // First, we need to determine the number of entries and the number of rows of the submatrix.
+            uint_fast64_t submatrixColumnCount = columnConstraint.getNumberOfSetBits();
+            
+            // Start by creating a temporary vector that stores for each index whose bit is set to true the number of
+            // bits that were set before that particular index.
+            std::vector<index_type> columnBitsSetBeforeIndex;
+            columnBitsSetBeforeIndex.reserve(columnCount);
+            
+            // Compute the information to fill this vector.
+            index_type lastIndex = 0;
+            index_type currentNumberOfSetBits = 0;
+            for (auto index : columnConstraint) {
+                while (lastIndex <= index) {
+                    columnBitsSetBeforeIndex.push_back(currentNumberOfSetBits);
+                    ++lastIndex;
+                }
+                ++currentNumberOfSetBits;
+            }
+            std::vector<index_type>* rowBitsSetBeforeIndex;
+            if (&rowGroupConstraint == &columnConstraint) {
+                rowBitsSetBeforeIndex = &columnBitsSetBeforeIndex;
+            } else {
+                rowBitsSetBeforeIndex = new std::vector<index_type>(rowCount);
+                
+                lastIndex = 0;
+                currentNumberOfSetBits = 0;
+                for (auto index : rowGroupConstraint) {
+                    while (lastIndex <= index) {
+                        rowBitsSetBeforeIndex->push_back(currentNumberOfSetBits);
+                        ++lastIndex;
+                    }
+                    ++currentNumberOfSetBits;
+                }
+            }
+            
+            // Then, we need to determine the number of entries and the number of rows of the submatrix.
             index_type subEntries = 0;
             index_type subRows = 0;
+            index_type rowGroupCount = 0;
             for (auto index : rowGroupConstraint) {
                 subRows += rowGroupIndices[index + 1] - rowGroupIndices[index];
                 for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
@@ -492,69 +560,55 @@ namespace storm {
                         if (columnConstraint.get(it->getColumn())) {
                             ++subEntries;
                             
-                            if (index == it->getColumn()) {
+                            if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) {
                                 foundDiagonalElement = true;
                             }
                         }
                     }
                     
                     // If requested, we need to reserve one entry more for inserting the diagonal zero entry.
-                    if (insertDiagonalEntries && !foundDiagonalElement) {
+                    if (insertDiagonalEntries && !foundDiagonalElement && rowGroupCount < submatrixColumnCount) {
                         ++subEntries;
                     }
                 }
+                ++rowGroupCount;
             }
             
             // Create and initialize resulting matrix.
-            SparseMatrixBuilder<ValueType> matrixBuilder(subRows, columnConstraint.getNumberOfSetBits(), subEntries, true, true);
-            
-            // Create a temporary vector that stores for each index whose bit is set to true the number of bits that
-            // were set before that particular index.
-            std::vector<index_type> bitsSetBeforeIndex;
-            bitsSetBeforeIndex.reserve(columnCount);
-            
-            // Compute the information to fill this vector.
-            index_type lastIndex = 0;
-            index_type currentNumberOfSetBits = 0;
-            
-            // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows are also
-            // taken.
-            storm::storage::BitVector columnBitCountConstraint = columnConstraint;
-            if (insertDiagonalEntries) {
-                columnBitCountConstraint |= rowGroupConstraint;
-            }
-            for (auto index : columnBitCountConstraint) {
-                while (lastIndex <= index) {
-                    bitsSetBeforeIndex.push_back(currentNumberOfSetBits);
-                    ++lastIndex;
-                }
-                ++currentNumberOfSetBits;
-            }
+            SparseMatrixBuilder<ValueType> matrixBuilder(subRows, submatrixColumnCount, subEntries, true, this->hasNontrivialRowGrouping());
             
             // Copy over selected entries.
+            rowGroupCount = 0;
             index_type rowCount = 0;
             for (auto index : rowGroupConstraint) {
-                matrixBuilder.newRowGroup(rowCount);
+                if (this->hasNontrivialRowGrouping()) {
+                    matrixBuilder.newRowGroup(rowCount);
+                }
                 for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
                     bool insertedDiagonalElement = false;
                     
                     for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) {
                         if (columnConstraint.get(it->getColumn())) {
-                            if (index == it->getColumn()) {
+                            if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) {
                                 insertedDiagonalElement = true;
-                            } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) {
-                                matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>());
+                            } else if (insertDiagonalEntries && !insertedDiagonalElement && columnBitsSetBeforeIndex[it->getColumn()] > (*rowBitsSetBeforeIndex)[index]) {
+                                matrixBuilder.addNextValue(rowCount, rowGroupCount, storm::utility::zero<ValueType>());
                                 insertedDiagonalElement = true;
                             }
-                            matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue());
+                            matrixBuilder.addNextValue(rowCount, columnBitsSetBeforeIndex[it->getColumn()], it->getValue());
                         }
                     }
-                    if (insertDiagonalEntries && !insertedDiagonalElement) {
-                        matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>());
+                    if (insertDiagonalEntries && !insertedDiagonalElement && rowGroupCount < submatrixColumnCount) {
+                        matrixBuilder.addNextValue(rowGroupCount, rowGroupCount, storm::utility::zero<ValueType>());
                     }
-                    
                     ++rowCount;
                 }
+                ++rowGroupCount;
+            }
+            
+            // If the constraints were not the same, we have allocated some additional memory we need to free now.
+            if (&rowGroupConstraint != &columnConstraint) {
+                delete rowBitsSetBeforeIndex;
             }
             
             return matrixBuilder.build();
@@ -615,12 +669,12 @@ namespace storm {
         SparseMatrix<ValueType> SparseMatrix<ValueType>::transpose(bool joinGroups, bool keepZeros) const {
             index_type rowCount = this->getColumnCount();
             index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount();
-			if (keepZeros) {
-				index_type entryCount = this->getEntryCount();
-			} else {
-				this->updateNonzeroEntryCount();
-				index_type entryCount = this->getNonzeroEntryCount();
-			}
+            if (keepZeros) {
+                index_type entryCount = this->getEntryCount();
+            } else {
+                this->updateNonzeroEntryCount();
+                index_type entryCount = this->getNonzeroEntryCount();
+            }
             
             std::vector<index_type> rowIndications(rowCount + 1);
             std::vector<MatrixEntry<index_type, ValueType>> columnsAndValues(entryCount);
@@ -628,7 +682,7 @@ namespace storm {
             // First, we need to count how many entries each column has.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         ++rowIndications[transition.getColumn() + 1];
                     }
                 }
@@ -647,7 +701,7 @@ namespace storm {
             // Now we are ready to actually fill in the values of the transposed matrix.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue());
                         nextIndices[transition.getColumn()]++;
                     }
@@ -659,11 +713,11 @@ namespace storm {
                 rowGroupIndices[i] = i;
             }
             
-            storm::storage::SparseMatrix<ValueType> transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices));
-                        
+            storm::storage::SparseMatrix<ValueType> transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), false);
+            
             return transposedMatrix;
         }
-                
+        
         template<typename ValueType>
         void SparseMatrix<ValueType>::convertToEquationSystem() {
             invertDiagonal();
@@ -673,22 +727,22 @@ namespace storm {
         template<typename ValueType>
         void SparseMatrix<ValueType>::invertDiagonal() {
             // Now iterate over all row groups and set the diagonal elements to the inverted value.
-			// If there is a row without the diagonal element, an exception is thrown.
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
+            // If there is a row without the diagonal element, an exception is thrown.
+            ValueType one = storm::utility::one<ValueType>();
+            ValueType zero = storm::utility::zero<ValueType>();
             bool foundDiagonalElement = false;
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
-						if (entry.getValue() == one) {
-							--this->nonzeroEntryCount;
-							entry.setValue(zero);
-						} else if (entry.getValue() == zero) {
-							++this->nonzeroEntryCount;
-							entry.setValue(one);
-						} else {
-							entry.setValue(one - entry.getValue());
-						}
+                        if (entry.getValue() == one) {
+                            --this->nonzeroEntryCount;
+                            entry.setValue(zero);
+                        } else if (entry.getValue() == zero) {
+                            ++this->nonzeroEntryCount;
+                            entry.setValue(one);
+                        } else {
+                            entry.setValue(one - entry.getValue());
+                        }
                         foundDiagonalElement = true;
                     }
                 }
@@ -718,7 +772,7 @@ namespace storm {
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
-						--this->nonzeroEntryCount;
+                        --this->nonzeroEntryCount;
                         entry.setValue(storm::utility::zero<ValueType>());
                     }
                 }
@@ -799,15 +853,28 @@ namespace storm {
             typename std::vector<ValueType>::iterator resultIterator = result.begin();
             typename std::vector<ValueType>::iterator resultIteratorEnd = result.end();
             
-            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
-                *resultIterator = storm::utility::zero<ValueType>();
-                
-                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
-                    *resultIterator += it->getValue() * vector[it->getColumn()];
+            // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable
+            // to store the intermediate result.
+            if (&vector == &result) {
+                for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                    ValueType tmpValue = storm::utility::zero<ValueType>();
+                    
+                    for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                        tmpValue += it->getValue() * vector[it->getColumn()];
+                    }
+                    *resultIterator = tmpValue;
+                }
+            } else {
+                for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                    *resultIterator = storm::utility::zero<ValueType>();
+                    
+                    for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                        *resultIterator += it->getValue() * vector[it->getColumn()];
+                    }
                 }
             }
         }
-
+        
 #ifdef STORM_HAVE_INTELTBB
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyWithVectorParallel(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const {
@@ -832,7 +899,53 @@ namespace storm {
                               });
         }
 #endif
-
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            const_iterator it = this->begin();
+            const_iterator ite;
+            std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
+            typename std::vector<ValueType>::const_iterator bIt = b.begin();
+            typename std::vector<ValueType>::const_iterator bIte = b.end();
+            typename std::vector<ValueType>::iterator resultIterator = x.begin();
+            typename std::vector<ValueType>::iterator resultIteratorEnd = x.end();
+            
+            // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable
+            // to store the intermediate result.
+            index_type currentRow = 0;
+            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++bIt) {
+                ValueType tmpValue = storm::utility::zero<ValueType>();
+                ValueType diagonalElement = storm::utility::zero<ValueType>();
+                
+                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                    if (it->getColumn() != currentRow) {
+                        tmpValue += it->getValue() * x[it->getColumn()];
+                    } else {
+                        diagonalElement += it->getValue();
+                    }
+                }
+                
+                *resultIterator = ((storm::utility::one<ValueType>() - omega) * *resultIterator) + (omega / diagonalElement) * (*bIt - tmpValue);
+                ++currentRow;
+            }
+        }
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const {
+            const_iterator it = this->begin();
+            const_iterator ite;
+            std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
+            std::vector<index_type>::const_iterator rowIteratorEnd = rowIndications.end();
+            
+            uint_fast64_t currentRow = 0;
+            for (; rowIterator != rowIteratorEnd - 1; ++rowIterator) {
+                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                    result[it->getColumn()] += it->getValue() * vector[currentRow];
+                }
+                ++currentRow;
+            }
+        }
+        
         template<typename ValueType>
         std::size_t SparseMatrix<ValueType>::getSizeInBytes() const {
             uint_fast64_t size = sizeof(*this);
@@ -906,6 +1019,11 @@ namespace storm {
             return this->columnsAndValues.begin() + this->rowIndications[rowCount];
         }
         
+        template<typename ValueType>
+        bool SparseMatrix<ValueType>::hasNontrivialRowGrouping() const {
+            return nontrivialRowGrouping;
+        }
+        
         template<typename ValueType>
         ValueType SparseMatrix<ValueType>::getRowSum(index_type row) const {
             ValueType sum = storm::utility::zero<ValueType>();
@@ -921,7 +1039,7 @@ namespace storm {
             if (this->getRowCount() != matrix.getRowCount()) return false;
             if (this->getColumnCount() != matrix.getColumnCount()) return false;
             if (this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false;
-
+            
             // Check the subset property for all rows individually.
             for (index_type row = 0; row < this->getRowCount(); ++row) {
                 for (const_iterator it1 = this->begin(row), ite1 = this->end(row), it2 = matrix.begin(row), ite2 = matrix.end(row); it1 != ite1; ++it1) {

src/storage/SparseMatrix.h

@@ -149,6 +149,14 @@ namespace storm {
              */
             SparseMatrixBuilder(index_type rows = 0, index_type columns = 0, index_type entries = 0, bool forceDimensions = true, bool hasCustomRowGrouping = false, index_type rowGroups = 0);
             
+            /*!
+             * Moves the contents of the given matrix into the matrix builder so that its contents can be modified again.
+             * This is, for example, useful if rows need to be added to the matrix.
+             *
+             * @param matrix The matrix that is to be made editable again.
+             */
+            SparseMatrixBuilder(SparseMatrix<ValueType>&& matrix);
+            
             /*!
              * Sets the matrix entry at the given row and column to the given value. After all entries have been added,
              * a call to finalize(false) is mandatory.
@@ -193,6 +201,20 @@ namespace storm {
              */
             SparseMatrix<value_type> build(index_type overriddenRowCount = 0, index_type overriddenColumnCount = 0, index_type overriddenRowGroupCount = 0);
             
+            /*!
+             * Retrieves the most recently used row.
+             * 
+             * @return The most recently used row.
+             */
+            index_type getLastRow() const;
+
+            /*!
+             * Retrieves the most recently used row.
+             *
+             * @return The most recently used row.
+             */
+            index_type getLastColumn() const;
+            
         private:
             // A flag indicating whether a row count was set upon construction.
             bool initialRowCountSet;
@@ -277,6 +299,7 @@ namespace storm {
             friend class storm::adapters::EigenAdapter;
             friend class storm::adapters::StormAdapter;
 			friend class storm::solver::TopologicalValueIterationMinMaxLinearEquationSolver<ValueType>;
+            friend class SparseMatrixBuilder<ValueType>;
             
             typedef uint_fast64_t index_type;
             typedef ValueType value_type;
@@ -388,6 +411,15 @@ namespace storm {
              */
             SparseMatrix(SparseMatrix<value_type> const& other);
             
+            /*!
+             * Constructs a sparse matrix by performing a deep-copy of the given matrix.
+             *
+             * @param other The matrix from which to copy the content.
+             * @param insertDiagonalElements If set to true, the copy will have all diagonal elements. If they did not
+             * exist in the original matrix, they are inserted and set to value zero.
+             */
+            SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements);
+            
             /*!
              * Constructs a sparse matrix by moving the contents of the given matrix to the newly created one.
              *
@@ -402,8 +434,9 @@ namespace storm {
              * @param rowIndications The row indications vector of the matrix to be constructed.
              * @param columnsAndValues The vector containing the columns and values of the entries in the matrix.
              * @param rowGroupIndices The vector representing the row groups in the matrix.
+             * @param hasNontrivialRowGrouping If set to true, this indicates that the row grouping is non-trivial.
              */
-            SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices);
+            SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool hasNontrivialRowGrouping);
             
             /*!
              * Constructs a sparse matrix by moving the given contents.
@@ -412,8 +445,9 @@ namespace storm {
              * @param rowIndications The row indications vector of the matrix to be constructed.
              * @param columnsAndValues The vector containing the columns and values of the entries in the matrix.
              * @param rowGroupIndices The vector representing the row groups in the matrix.
+             * @param hasNontrivialRowGrouping If set to true, this indicates that the row grouping is non-trivial.
              */
-            SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, value_type>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices);
+            SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, value_type>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices, bool hasNontrivialRowGrouping);
 
             /*!
              * Assigns the contents of the given matrix to the current one by deep-copying its contents.
@@ -646,6 +680,26 @@ namespace storm {
              * @return The product of the matrix and the given vector as the content of the given result vector.
              */
             void multiplyWithVector(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
+
+            /*!
+             * Multiplies the vector to the matrix from the left and writes the result to the given result vector.
+             *
+             * @param vector The vector with which the matrix is to be multiplied. This vector is interpreted as being
+             * a row vector.
+             * @param result The vector that is supposed to hold the result of the multiplication after the operation.
+             * @return The product of the matrix and the given vector as the content of the given result vector. The 
+             * result is to be interpreted as a row vector.
+             */
+            void multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
+            
+            /*!
+             * Performs one step of the successive over-relaxation technique.
+             *
+             * @param omega The Omega parameter for SOR.
+             * @param x The current solution vector. The result will be written to the very same vector.
+             * @param b The 'right-hand side' of the problem.
+             */
+            void performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             
             /*!
              * Multiplies the matrix with the given vector in a sequential way and writes the result to the given result
@@ -798,6 +852,13 @@ namespace storm {
              * @return An iterator that points past the end of the last row of the matrix.
              */
             iterator end();
+            
+            /*!
+             * Retrieves whether the matrix has a (possibly) non-trivial row grouping.
+             *
+             * @return True iff the matrix has a (possibly) non-trivial row grouping.
+             */
+            bool hasNontrivialRowGrouping() const;
 
 			/*!
 			* Returns a copy of the matrix with the chosen internal data type
@@ -813,7 +874,7 @@ namespace storm {
 					new_columnsAndValues.at(i) = MatrixEntry<SparseMatrix::index_type, NewValueType>(columnsAndValues.at(i).getColumn(), static_cast<NewValueType>(columnsAndValues.at(i).getValue()));
 				}
 
-				return SparseMatrix<NewValueType>(columnCount, std::move(new_rowIndications), std::move(new_columnsAndValues), std::move(new_rowGroupIndices));
+				return SparseMatrix<NewValueType>(columnCount, std::move(new_rowIndications), std::move(new_columnsAndValues), std::move(new_rowGroupIndices), nontrivialRowGrouping);
 			}
             
         private:
@@ -852,6 +913,10 @@ namespace storm {
             // entry is not included anymore.
             std::vector<index_type> rowIndications;
             
+            // A flag that indicates whether the matrix has a non-trivial row-grouping, i.e. (possibly) more than one
+            // row per row group.
+            bool nontrivialRowGrouping;
+            
             // A vector indicating the row groups of the matrix.
             std::vector<index_type> rowGroupIndices;
         };

src/storage/dd/CuddAdd.cpp

@@ -557,7 +557,7 @@ namespace storm {
             rowIndications[0] = 0;
             
             // Construct matrix and return result.
-            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(trivialRowGroupIndices));
+            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(trivialRowGroupIndices), false);
         }
         
         storm::storage::SparseMatrix<double> Add<DdType::CUDD>::toMatrix(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
@@ -680,7 +680,7 @@ namespace storm {
             }
             rowIndications[0] = 0;
             
-            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices));
+            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), true);
         }
         
         std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> Add<DdType::CUDD>::toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
@@ -804,7 +804,7 @@ namespace storm {
             }
             rowIndications[0] = 0;
             
-            return std::make_pair(storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)), std::move(explicitVector));
+            return std::make_pair(storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), true), std::move(explicitVector));
         }
         
         template<typename ValueType>

src/storage/dd/CuddAdd.h

@@ -721,7 +721,7 @@ namespace storm {
             void toMatrixRec(DdNode const* dd, std::vector<uint_fast64_t>& rowIndications, std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>>& columnsAndValues, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, Odd<DdType::CUDD> const& columnOdd, uint_fast64_t currentRowLevel, uint_fast64_t currentColumnLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, uint_fast64_t currentColumnOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices, std::vector<uint_fast64_t> const& ddColumnVariableIndices, bool generateValues = true) const;
             
             /*!
-             * Helper function to convert the DD into a (sparse) vector.
+             * Helper function to convert the DD into an (explicit) vector.
              *
              * @param dd The DD to convert.
              * @param result The vector that will hold the values upon successful completion.
@@ -772,7 +772,7 @@ namespace storm {
             void addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const;
             
             /*!
-             * Performs a recursive step to perform the given function between the given DD-based vector to the given
+             * Performs a recursive step to perform the given function between the given DD-based vector and the given
              * explicit vector.
              *
              * @param dd The DD to add to the explicit vector.

src/storage/dd/CuddBdd.cpp

@@ -337,6 +337,41 @@ namespace storm {
             }
         }
         
+        storm::storage::BitVector Bdd<DdType::CUDD>::toVector(storm::dd::Odd<DdType::CUDD> const& rowOdd) const {
+            std::vector<uint_fast64_t> ddVariableIndices = this->getSortedVariableIndices();
+            storm::storage::BitVector result(rowOdd.getTotalOffset());
+            this->toVectorRec(this->getCuddDdNode(), this->getDdManager()->getCuddManager(), result, rowOdd, Cudd_IsComplement(this->getCuddDdNode()), 0, ddVariableIndices.size(), 0, ddVariableIndices);
+            return result;
+        }
+        
+        void Bdd<DdType::CUDD>::toVectorRec(DdNode const* dd, Cudd const& manager, storm::storage::BitVector& result, Odd<DdType::CUDD> const& rowOdd, bool complement, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const {
+            // If there are no more values to select, we can directly return.
+            if (dd == Cudd_ReadLogicZero(manager.getManager()) && !complement) {
+                return;
+            } else if (dd == Cudd_ReadOne(manager.getManager()) && complement) {
+                return;
+            }
+            
+            // If we are at the maximal level, the value to be set is stored as a constant in the DD.
+            if (currentRowLevel == maxLevel) {
+                result.set(currentRowOffset, true);
+            } else if (ddRowVariableIndices[currentRowLevel] < dd->index) {
+                toVectorRec(dd, manager, result, rowOdd.getElseSuccessor(), complement, currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices);
+                toVectorRec(dd, manager, result, rowOdd.getThenSuccessor(), complement, currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices);
+            } else {
+                // Otherwise, we compute the ODDs for both the then- and else successors.
+                DdNode* elseDdNode = Cudd_E(dd);
+                DdNode* thenDdNode = Cudd_T(dd);
+                
+                // Determine whether we have to evaluate the successors as if they were complemented.
+                bool elseComplemented = Cudd_IsComplement(elseDdNode) ^ complement;
+                bool thenComplemented = Cudd_IsComplement(thenDdNode) ^ complement;
+                
+                toVectorRec(Cudd_Regular(elseDdNode), manager, result, rowOdd.getElseSuccessor(), elseComplemented, currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices);
+                toVectorRec(Cudd_Regular(thenDdNode), manager, result, rowOdd.getThenSuccessor(), thenComplemented, currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices);
+            }
+        }
+        
         std::ostream& operator<<(std::ostream& out, const Bdd<DdType::CUDD>& bdd) {
             bdd.exportToDot();
             return out;

src/storage/dd/CuddBdd.h

@@ -263,6 +263,15 @@ namespace storm {
              */
             Add<DdType::CUDD> toAdd() const;
             
+            /*!
+             * Converts the BDD to a bit vector. The given offset-labeled DD is used to determine the correct row of
+             * each entry.
+             *
+             * @param rowOdd The ODD used for determining the correct row.
+             * @return The bit vector that is represented by this BDD.
+             */
+            storm::storage::BitVector toVector(storm::dd::Odd<DdType::CUDD> const& rowOdd) const;
+            
         private:
             /*!
              * Retrieves the CUDD BDD object associated with this DD.
@@ -315,6 +324,21 @@ namespace storm {
             template<typename ValueType>
             static DdNode* fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::function<bool (ValueType const&)> const& filter);
             
+            /*!
+             * Helper function to convert the DD into a bit vector.
+             *
+             * @param dd The DD to convert.
+             * @param manager The Cudd manager responsible for the DDs.
+             * @param result The vector that will hold the values upon successful completion.
+             * @param rowOdd The ODD used for the row translation.
+             * @param complement A flag indicating whether the result is to be interpreted as a complement.
+             * @param currentRowLevel The currently considered row level in the DD.
+             * @param maxLevel The number of levels that need to be considered.
+             * @param currentRowOffset The current row offset.
+             * @param ddRowVariableIndices The (sorted) indices of all DD row variables that need to be considered.
+             */
+            void toVectorRec(DdNode const* dd, Cudd const& manager, storm::storage::BitVector& result, Odd<DdType::CUDD> const& rowOdd, bool complement, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const;
+            
             // The BDD created by CUDD.
             BDD cuddBdd;
         };

src/storage/expressions/ExprtkExpressionEvaluator.cpp

@@ -49,7 +49,7 @@ namespace storm {
             CompiledExpressionType& compiledExpression = result.first->second;
             compiledExpression.register_symbol_table(symbolTable);
             bool parsingOk = parser.compile(ToExprtkStringVisitor().toString(expression), compiledExpression);
-            STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk.");
+            STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk: " << *expression);
             return compiledExpression;
         }
         

src/utility/cli.cpp

@@ -34,11 +34,11 @@ namespace storm {
             }
             
             void printHeader(const int argc, const char* argv[]) {
-                std::cout << "StoRM" << std::endl;
+                 std::cout << "StoRM" << std::endl;
                 std::cout << "--------" << std::endl << std::endl;
                 
                 
-                //				std::cout << storm::utility::StormVersion::longVersionString() << std::endl;
+                std::cout << storm::utility::StormVersion::longVersionString() << std::endl;
 #ifdef STORM_HAVE_INTELTBB
                 std::cout << "Linked with Intel Threading Building Blocks v" << TBB_VERSION_MAJOR << "." << TBB_VERSION_MINOR << " (Interface version " << TBB_INTERFACE_VERSION << ")." << std::endl;
 #endif
@@ -113,6 +113,7 @@ namespace storm {
                 std::cout << "Command line arguments: " << commandStream.str() << std::endl;
                 std::cout << "Current working directory: " << getCurrentWorkingDirectory() << std::endl << std::endl;
             }
+            }
             
             
             void printUsage() {

src/utility/cli.h

@@ -251,6 +251,22 @@ namespace storm {
             }
             
 #ifdef STORM_HAVE_CARL
+            void exportParametricResultToFile(storm::RationalFunction const& result, storm::models::sparse::Dtmc<storm::RationalFunction>::ConstraintCollector const& constraintCollector, std::string const& path) {
+                std::ofstream filestream;
+                filestream.open(path);
+                // TODO: add checks.
+                filestream << "!Parameters: ";
+                std::set<storm::Variable> vars = result.gatherVariables();
+                std::copy(vars.begin(), vars.end(), std::ostream_iterator<storm::Variable>(filestream, ", "));
+                filestream << std::endl;
+                filestream << "!Result: " << result << std::endl;
+                filestream << "!Well-formed Constraints: " << std::endl;
+                std::copy(constraintCollector.getWellformedConstraints().begin(), constraintCollector.getWellformedConstraints().end(), std::ostream_iterator<carl::Constraint<storm::RationalFunction>>(filestream, "\n"));
+                filestream << "!Graph-preserving Constraints: " << std::endl;
+                std::copy(constraintCollector.getGraphPreservingConstraints().begin(), constraintCollector.getGraphPreservingConstraints().end(), std::ostream_iterator<carl::Constraint<storm::RationalFunction>>(filestream, "\n"));
+                filestream.close();
+            }
+            
             template<>
             inline void verifySparseModel(boost::optional<storm::prism::Program> const& program, std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>> model, std::shared_ptr<storm::logic::Formula> formula) {
                 storm::settings::modules::GeneralSettings const& settings = storm::settings::generalSettings();
@@ -276,6 +292,11 @@ namespace storm {
                 } else {
                     std::cout << " skipped, because the modelling formalism is currently unsupported." << std::endl;
                 }
+                
+                storm::settings::modules::ParametricSettings const& parametricSettings = storm::settings::parametricSettings();
+                if (parametricSettings.exportResultToFile()) {
+                    exportParametricResultToFile(result->asExplicitQuantitativeCheckResult<storm::RationalFunction>()[*dtmc->getInitialStates().begin()], storm::models::sparse::Dtmc<storm::RationalFunction>::ConstraintCollector(*dtmc), parametricSettings.exportResultPath());
+                }
             }
 #endif
             
@@ -462,4 +483,4 @@ namespace storm {
     }
 }
 
-#endif
+#endif

src/utility/graph.h

@@ -245,7 +245,7 @@ namespace storm {
              * with probability 0 and the second stores all indices of states with probability 1.
              */
             template <typename T>
-            static std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01(storm::storage::SparseMatrix<T> backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+            static std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01(storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
                 std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
                 result.first = performProbGreater0(backwardTransitions, phiStates, psiStates);
                 result.second = performProb1(backwardTransitions, phiStates, psiStates, result.first);

src/utility/solver.cpp

@@ -2,6 +2,8 @@
 
 #include "src/settings/SettingsManager.h"
 
+#include "src/solver/SymbolicGameSolver.h"
+
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
 
@@ -15,6 +17,16 @@
 namespace storm {
     namespace utility {
         namespace solver {
+            template<storm::dd::DdType Type, typename ValueType>
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> SymbolicLinearEquationSolverFactory<Type, ValueType>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
+                return std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicLinearEquationSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs));
+            }
+            
+            template<storm::dd::DdType Type>
+            std::unique_ptr<storm::solver::SymbolicGameSolver<Type>> SymbolicGameSolverFactory<Type>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const {
+                return std::unique_ptr<storm::solver::SymbolicGameSolver<Type>>(new storm::solver::SymbolicGameSolver<Type>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables));
+            }
+            
             template<typename ValueType>
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
                 storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver();
@@ -29,9 +41,28 @@ namespace storm {
                 return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
             }
             
+            template<typename ValueType>
+            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory() {
+                switch (storm::settings::nativeEquationSolverSettings().getLinearEquationSystemMethod()) {
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::Jacobi;
+                    break;
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::GaussSeidel;
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::SOR;
+                }
+                omega = storm::settings::nativeEquationSolverSettings().getOmega();
+            }
+            
+            template<typename ValueType>
+            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega) : method(method), omega(omega) {
+                // Intentionally left empty.
+            }
+            
             template<typename ValueType>
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
+                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix, method));
             }
             
             template<typename ValueType>
@@ -79,6 +110,8 @@ namespace storm {
                 return factory->create(name);
             }
             
+            template class SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>;
+            template class SymbolicGameSolverFactory<storm::dd::DdType::CUDD>;
             template class LinearEquationSolverFactory<double>;
             template class GmmxxLinearEquationSolverFactory<double>;
             template class NativeLinearEquationSolverFactory<double>;

src/utility/solver.h

@@ -1,15 +1,34 @@
 #ifndef STORM_UTILITY_SOLVER_H_
 #define STORM_UTILITY_SOLVER_H_
 
+#include "src/solver/SymbolicGameSolver.h"
+
+#include "src/solver/SymbolicLinearEquationSolver.h"
 #include "src/solver/LinearEquationSolver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/MinMaxLinearEquationSolver.h"
 #include "src/solver/LpSolver.h"
 
+#include "src/storage/dd/DdType.h"
+#include "src/settings/modules/NativeEquationSolverSettings.h"
+
 #include "src/exceptions/InvalidSettingsException.h"
 
 namespace storm {
     namespace utility {
         namespace solver {
+            template<storm::dd::DdType Type, typename ValueType>
+            class SymbolicLinearEquationSolverFactory {
+            public:
+                virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
+            };
+            
+            template<storm::dd::DdType Type>
+            class SymbolicGameSolverFactory {
+            public:
+                virtual std::unique_ptr<storm::solver::SymbolicGameSolver<Type>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const;
+            };
+            
             template<typename ValueType>
             class LinearEquationSolverFactory {
             public:
@@ -25,7 +44,14 @@ namespace storm {
             template<typename ValueType>
             class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
             public:
+                NativeLinearEquationSolverFactory();
+                NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega);
+                
                 virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+                
+            private:
+                typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method;
+                ValueType omega;
             };
             
             template<typename ValueType>

src/utility/vector.h

@@ -376,7 +376,9 @@ namespace storm {
 					if (val2 == 0) {
 						return (std::abs(val1) <= precision);
 					}
-                    if (std::abs((val1 - val2)/val2) > precision) return false;
+                    if (std::abs((val1 - val2)/val2) > precision) {
+                        return false;
+                    }
                 } else {
                     if (std::abs(val1 - val2) > precision) return false;
                 }

test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp

@@ -86,6 +86,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) {
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision());
+
+    formula = formulaParser.parseFromString("LRA=? [\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult8 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
@@ -149,6 +156,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+
+    formula = formulaParser.parseFromString("LRA=? [\"fail_sensors\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.93458866427696596, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxCtmcCslModelCheckerTest, Polling) {
@@ -176,6 +190,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Polling) {
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult1 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(1, quantitativeCheckResult1[initialState], storm::settings::generalSettings().getPrecision());
+
+    formula = formulaParser.parseFromString("LRA=?[\"target\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult2 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxCtmcCslModelCheckerTest, Fms) {
@@ -252,4 +273,11 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) {
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(262.85103661561755, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("LRA=? [\"first_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.9100373532, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision());
 }

test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp

@@ -127,3 +127,148 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) {
 
 	EXPECT_NEAR(1.044879046, quantitativeResult3[0], storm::settings::gmmxxEquationSolverSettings().getPrecision());
 }
+
+TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2);
+        matrixBuilder.addNextValue(0, 1, 1.);
+        matrixBuilder.addNextValue(1, 0, 1.);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+        matrixBuilder.addNextValue(0, 0, .5);
+        matrixBuilder.addNextValue(0, 1, .5);
+        matrixBuilder.addNextValue(1, 0, .5);
+        matrixBuilder.addNextValue(1, 1, .5);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 2, 1);
+        matrixBuilder.addNextValue(2, 0, 1);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(3);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 2);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+}
+
+TEST(GmmxxDtmcPrctlModelCheckerTest, LRA) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> mdp;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(15, 15, 20, true);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 4, 0.7);
+        matrixBuilder.addNextValue(1, 6, 0.3);
+        matrixBuilder.addNextValue(2, 0, 1);
+        
+        matrixBuilder.addNextValue(3, 5, 0.8);
+        matrixBuilder.addNextValue(3, 9, 0.2);
+        matrixBuilder.addNextValue(4, 3, 1);
+        matrixBuilder.addNextValue(5, 3, 1);
+        
+        matrixBuilder.addNextValue(6, 7, 1);
+        matrixBuilder.addNextValue(7, 8, 1);
+        matrixBuilder.addNextValue(8, 6, 1);
+        
+        matrixBuilder.addNextValue(9, 10, 1);
+        matrixBuilder.addNextValue(10, 9, 1);
+        matrixBuilder.addNextValue(11, 9, 1);
+        
+        matrixBuilder.addNextValue(12, 5, 0.4);
+        matrixBuilder.addNextValue(12, 8, 0.3);
+        matrixBuilder.addNextValue(12, 11, 0.3);
+        
+        matrixBuilder.addNextValue(13, 7, 0.7);
+        matrixBuilder.addNextValue(13, 12, 0.3);
+        
+        matrixBuilder.addNextValue(14, 12, 1);
+        
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(15);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        ap.addLabelToState("a", 4);
+        ap.addLabelToState("a", 5);
+        ap.addLabelToState("a", 7);
+        ap.addLabelToState("a", 11);
+        ap.addLabelToState("a", 13);
+        ap.addLabelToState("a", 14);
+        
+        mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.3/3., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+}

test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp

@@ -102,6 +102,14 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Cluster) {
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("LRA=? [\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult8 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8.getMax(), storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded) {
@@ -173,6 +181,14 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded) {
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("LRA=? [\"fail_sensors\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.934586179, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.934586179, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling) {
@@ -201,6 +217,15 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling) {
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     EXPECT_NEAR(1, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(1, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("LRA=? [\"target\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+
 }
 
 TEST(GmmxxHybridCtmcCslModelCheckerTest, Fms) {
@@ -288,4 +313,12 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Tandem) {
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("LRA=? [\"first_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.9100373532, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9100373532, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision());
 }

test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp

@@ -1,6 +1,7 @@
 #include "gtest/gtest.h"
 #include "storm-config.h"
 
+#include "src/settings/SettingMemento.h"
 #include "src/logic/Formulas.h"
 #include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
@@ -128,7 +129,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, SynchronousLeader) {
     EXPECT_NEAR(1.0448979589010925, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
-TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
+TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) {
 	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
 	std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
 
@@ -144,7 +145,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -169,7 +170,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -194,7 +195,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -255,7 +256,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA) {
 
 		mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);

test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp

@@ -0,0 +1,172 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/settings/SettingsManager.h"
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, Die) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/die.pm");
+    
+    // Build the die model with its reward model.
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
+    options.buildRewards = true;
+    options.rewardModelName = "coin_flips";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(13, model->getNumberOfStates());
+    EXPECT_EQ(20, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("one");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto done = std::make_shared<storm::logic::AtomicLabelFormula>("done");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(done);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(3.6666622161865234, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(3.6666622161865234, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, Crowds) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/crowds-5-5.pm");
+    
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program);
+    EXPECT_EQ(8607, model->getNumberOfStates());
+    EXPECT_EQ(15113, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observe0Greater1");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.33288236360191303, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.33288236360191303, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeIGreater1");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.15222081144084315, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.15222081144084315, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeOnlyTrueSender");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.3215392962289586, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.3215392962289586, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, SynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader-3-5.pm");
+    
+    // Build the die model with its reward model.
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
+    options.buildRewards = true;
+    options.rewardModelName = "num_rounds";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(273, model->getNumberOfStates());
+    EXPECT_EQ(397, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 20);
+    
+    result = checker.check(*boundedUntilFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+

test/functional/solver/FullySymbolicGameSolverTest.cpp

@@ -0,0 +1,64 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/storage/dd/CuddDdManager.h"
+
+#include "src/utility/solver.h"
+
+TEST(FullySymbolicGameSolverTest, Solve) {
+    // Create some variables.
+    std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::CUDD>> manager(new storm::dd::DdManager<storm::dd::DdType::CUDD>());
+    std::pair<storm::expressions::Variable, storm::expressions::Variable> state = manager->addMetaVariable("x", 1, 4);
+    std::pair<storm::expressions::Variable, storm::expressions::Variable> pl1 = manager->addMetaVariable("a", 0, 1);
+    std::pair<storm::expressions::Variable, storm::expressions::Variable> pl2 = manager->addMetaVariable("b", 0, 1);
+
+    storm::dd::Bdd<storm::dd::DdType::CUDD> allRows = manager->getBddZero();
+    std::set<storm::expressions::Variable> rowMetaVariables({state.first});
+    std::set<storm::expressions::Variable> columnMetaVariables({state.second});
+    std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> rowColumnMetaVariablePairs = {state};
+    std::set<storm::expressions::Variable> player1Variables({pl1.first});
+    std::set<storm::expressions::Variable> player2Variables({pl2.first});
+
+    // Construct simple game.
+    storm::dd::Add<storm::dd::DdType::CUDD> matrix = manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 2).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.6);
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 1).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.4);
+
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 2).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(0.2);
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 3).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(0.8);
+
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 3).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.5);
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 4).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.5);
+    
+    matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 1).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(1);
+    
+    std::unique_ptr<storm::utility::solver::SymbolicGameSolverFactory<storm::dd::DdType::CUDD>> solverFactory(new storm::utility::solver::SymbolicGameSolverFactory<storm::dd::DdType::CUDD>());
+    std::unique_ptr<storm::solver::SymbolicGameSolver<storm::dd::DdType::CUDD>> solver = solverFactory->create(matrix, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables,player2Variables);
+    
+    // Create solution and target state vector.
+    storm::dd::Add<storm::dd::DdType::CUDD> x = manager->getAddZero();
+    storm::dd::Add<storm::dd::DdType::CUDD> b = manager->getEncoding(state.first, 2).toAdd() + manager->getEncoding(state.first, 4).toAdd();
+    
+    // Now solve the game with different strategies for the players.
+    storm::dd::Add<storm::dd::DdType::CUDD> result = solver->solveGame(true, true, x, b);
+    result *= manager->getEncoding(state.first, 1).toAdd();
+    result = result.sumAbstract({state.first});
+    EXPECT_NEAR(0, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    x = manager->getAddZero();
+    result = solver->solveGame(true, false, x, b);
+    result *= manager->getEncoding(state.first, 1).toAdd();
+    result = result.sumAbstract({state.first});
+    EXPECT_NEAR(0.5, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    x = manager->getAddZero();
+    result = solver->solveGame(false, true, x, b);
+    result *= manager->getEncoding(state.first, 1).toAdd();
+    result = result.sumAbstract({state.first});
+    EXPECT_NEAR(0.2, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision());
+
+    x = manager->getAddZero();
+    result = solver->solveGame(false, false, x, b);
+    result *= manager->getEncoding(state.first, 1).toAdd();
+    result = result.sumAbstract({state.first});
+    EXPECT_NEAR(0.99999892625817599, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}

test/functional/storage/SparseMatrixTest.cpp

@@ -154,8 +154,8 @@ TEST(SparseMatrix, CreationWithMovingContents) {
     columnsAndValues.emplace_back(1, 0.7);
     columnsAndValues.emplace_back(3, 0.2);
     
-    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}));
-    storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3});
+    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}, false));
+    storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}, false);
     ASSERT_EQ(3, matrix.getRowCount());
     ASSERT_EQ(4, matrix.getColumnCount());
     ASSERT_EQ(5, matrix.getEntryCount());
@@ -328,7 +328,7 @@ TEST(SparseMatrix, Submatrix) {
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 3, 0.3));
     storm::storage::SparseMatrix<double> matrix;
     ASSERT_NO_THROW(matrix = matrixBuilder.build());
-    
+
     std::vector<uint_fast64_t> rowGroupIndices = {0, 1, 2, 4, 5};
     
     storm::storage::BitVector rowGroupConstraint(4);