diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 88e0a396f..05d262b65 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -33,6 +33,7 @@ file(GLOB STORM_SETTINGS_FILES ${PROJECT_SOURCE_DIR}/src/settings/*.h ${PROJECT_
 file(GLOB STORM_SETTINGS_MODULES_FILES ${PROJECT_SOURCE_DIR}/src/settings/modules/*.h ${PROJECT_SOURCE_DIR}/src/settings/modules/*.cpp)
 file(GLOB_RECURSE STORM_SOLVER_FILES ${PROJECT_SOURCE_DIR}/src/solver/*.h ${PROJECT_SOURCE_DIR}/src/solver/*.cpp)
 file(GLOB STORM_STORAGE_FILES ${PROJECT_SOURCE_DIR}/src/storage/*.h ${PROJECT_SOURCE_DIR}/src/storage/*.cpp)
+file(GLOB STORM_STORAGE_BISIMULATION_FILES ${PROJECT_SOURCE_DIR}/src/storage/bisimulation/*.h ${PROJECT_SOURCE_DIR}/src/storage/bisimulation/*.cpp)
 file(GLOB_RECURSE STORM_STORAGE_DD_FILES ${PROJECT_SOURCE_DIR}/src/storage/dd/*.h ${PROJECT_SOURCE_DIR}/src/storage/dd/*.cpp)
 file(GLOB_RECURSE STORM_STORAGE_EXPRESSIONS_FILES ${PROJECT_SOURCE_DIR}/src/storage/expressions/*.h ${PROJECT_SOURCE_DIR}/src/storage/expressions/*.cpp)
 file(GLOB_RECURSE STORM_STORAGE_PRISM_FILES ${PROJECT_SOURCE_DIR}/src/storage/prism/*.h ${PROJECT_SOURCE_DIR}/src/storage/prism/*.cpp)
@@ -77,6 +78,7 @@ source_group(settings FILES ${STORM_SETTINGS_FILES})
 source_group(settings\\modules FILES ${STORM_SETTINGS_MODULES_FILES})
 source_group(solver FILES ${STORM_SOLVER_FILES})
 source_group(storage FILES ${STORM_STORAGE_FILES})
+source_group(storage\\bisimulation FILES ${STORM_STORAGE_BISIMULATION_FILES})
 source_group(storage\\dd FILES ${STORM_STORAGE_DD_FILES})
 source_group(storage\\expressions FILES ${STORM_STORAGE_EXPRESSIONS_FILES})
 source_group(storage\\prism FILES ${STORM_STORAGE_PRISM_FILES})
diff --git a/src/adapters/CarlAdapter.h b/src/adapters/CarlAdapter.h
index 99e01e575..2ed172cde 100644
--- a/src/adapters/CarlAdapter.h
+++ b/src/adapters/CarlAdapter.h
@@ -42,7 +42,7 @@ namespace carl {
 }
 
 namespace storm {
-    typedef mpq_class RationalNumber;
+    typedef cln::cl_RA RationalNumber;
     typedef carl::Variable Variable;
     typedef carl::MultivariatePolynomial<RationalNumber> RawPolynomial;
     typedef carl::FactorizedPolynomial<RawPolynomial> Polynomial;
diff --git a/src/builder/DdPrismModelBuilder.cpp b/src/builder/DdPrismModelBuilder.cpp
index 4f6ab3871..0a20bf38b 100644
--- a/src/builder/DdPrismModelBuilder.cpp
+++ b/src/builder/DdPrismModelBuilder.cpp
@@ -1020,7 +1020,7 @@ namespace storm {
             
             preparedProgram = preparedProgram.substituteConstants();
             
-            STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl);
+            STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl);
             
             // Start by initializing the structure used for storing all information needed during the model generation.
             // In particular, this creates the meta variables used to encode the model.
diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp
index 5b62364b4..4f6117c9c 100644
--- a/src/builder/ExplicitPrismModelBuilder.cpp
+++ b/src/builder/ExplicitPrismModelBuilder.cpp
@@ -295,7 +295,7 @@ namespace storm {
             // Now that the program is fixed, we we need to substitute all constants with their concrete value.
             preparedProgram = preparedProgram.substituteConstants();
             
-            STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl);
+            STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl);
                 
             // Select the appropriate reward models (after the constants have been substituted).
             std::vector<std::reference_wrapper<storm::prism::RewardModel const>> selectedRewardModels;
diff --git a/src/builder/ExplicitPrismModelBuilder.h b/src/builder/ExplicitPrismModelBuilder.h
index 8680c0d7b..dfab62136 100644
--- a/src/builder/ExplicitPrismModelBuilder.h
+++ b/src/builder/ExplicitPrismModelBuilder.h
@@ -72,8 +72,6 @@ namespace storm {
                 std::string stateInfo(uint_fast64_t state) const override {
                     return valuations[state].toString();
                 }
-
-
             };
             
             // A structure storing information about the used variables of the program.
diff --git a/src/cli/entrypoints.h b/src/cli/entrypoints.h
index c17b2fda9..bcbf8cdfa 100644
--- a/src/cli/entrypoints.h
+++ b/src/cli/entrypoints.h
@@ -3,11 +3,9 @@
 
 #include "src/utility/storm.h"
 
-
 namespace storm {
     namespace cli {
 
-
         template<typename ValueType>
         void verifySparseModel(std::shared_ptr<storm::models::sparse::Model<ValueType>> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
             for (auto const& formula : formulas) {
@@ -85,15 +83,40 @@ namespace storm {
                 }
             }
         }
-
-
+        
+#define BRANCH_ON_MODELTYPE(result, model, value_type, dd_type, function, ...) \
+    if (model->isSymbolicModel()) { \
+        if (model->isOfType(storm::models::ModelType::Dtmc)) { \
+            result = function<storm::models::symbolic::Dtmc<dd_type>>(model->as<storm::models::symbolic::Dtmc<dd_type>>(), __VA_ARGS__); \
+        } else if (model->isOfType(storm::models::ModelType::Ctmc)) { \
+            result = function<storm::models::symbolic::Ctmc<dd_type>>(model->as<storm::models::symbolic::Ctmc<dd_type>>(), __VA_ARGS__); \
+        } else if (model->isOfType(storm::models::ModelType::Mdp)) { \
+            result = function<storm::models::symbolic::Mdp<dd_type>>(model->as<storm::models::symbolic::Mdp<dd_type>>(), __VA_ARGS__); \
+        } else { \
+            STORM_LOG_ASSERT(false, "Unknown model type."); \
+        } \
+    } else { \
+        STORM_LOG_ASSERT(model->isSparseModel(), "Unknown model type."); \
+        if (model->isOfType(storm::models::ModelType::Dtmc)) { \
+            result = function<storm::models::sparse::Dtmc<value_type>>(model->as<storm::models::sparse::Dtmc<value_type>>(), __VA_ARGS__); \
+        } else if (model->isOfType(storm::models::ModelType::Ctmc)) { \
+            result = function<storm::models::sparse::Ctmc<value_type>>(model->as<storm::models::sparse::Ctmc<value_type>>(), __VA_ARGS__); \
+        } else if (model->isOfType(storm::models::ModelType::Mdp)) { \
+            result = function<storm::models::sparse::Mdp<value_type>>(model->as<storm::models::sparse::Mdp<value_type>>(), __VA_ARGS__); \
+        } else if (model->isOfType(storm::models::ModelType::MarkovAutomaton)) { \
+            result = function<storm::models::sparse::MarkovAutomaton<value_type>>(model->as<storm::models::sparse::MarkovAutomaton<value_type>>(), __VA_ARGS__); \
+        } else { \
+            STORM_LOG_ASSERT(false, "Unknown model type."); \
+        } \
+    }
+    
         template<typename ValueType>
         void buildAndCheckSymbolicModel(storm::prism::Program const& program, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
             std::shared_ptr<storm::models::ModelBase> model = buildSymbolicModel<ValueType>(program, formulas);
             STORM_LOG_THROW(model != nullptr, storm::exceptions::InvalidStateException, "Model could not be constructed for an unknown reason.");
 
             // Preprocess the model if needed.
-            model = preprocessModel<ValueType>(model, formulas);
+            BRANCH_ON_MODELTYPE(model, model, ValueType, storm::dd::DdType::CUDD, preprocessModel, formulas);
 
             // Print some information about the model.
             model->printModelInformationToStream(std::cout);
@@ -127,8 +150,9 @@ namespace storm {
 
             STORM_LOG_THROW(settings.isExplicitSet(), storm::exceptions::InvalidStateException, "Unable to build explicit model without model files.");
             std::shared_ptr<storm::models::ModelBase> model = buildExplicitModel<ValueType>(settings.getTransitionFilename(), settings.getLabelingFilename(), settings.isStateRewardsSet() ? settings.getStateRewardsFilename() : boost::optional<std::string>(), settings.isTransitionRewardsSet() ? settings.getTransitionRewardsFilename() : boost::optional<std::string>(), settings.isChoiceLabelingSet() ? settings.getChoiceLabelingFilename() : boost::optional<std::string>());
+            
             // Preprocess the model if needed.
-            model = preprocessModel<ValueType>(model, formulas);
+            BRANCH_ON_MODELTYPE(model, model, ValueType, storm::dd::DdType::CUDD, preprocessModel, formulas);
 
             // Print some information about the model.
             model->printModelInformationToStream(std::cout);
diff --git a/src/models/ModelBase.cpp b/src/models/ModelBase.cpp
index 1eff87d18..b79d7a81f 100644
--- a/src/models/ModelBase.cpp
+++ b/src/models/ModelBase.cpp
@@ -13,5 +13,9 @@ namespace storm {
         bool ModelBase::isSymbolicModel() const {
             return false;
         }
+        
+        bool ModelBase::isOfType(storm::models::ModelType const& modelType) const {
+            return this->getType() == modelType;
+        }
     }
 }
\ No newline at end of file
diff --git a/src/models/ModelBase.h b/src/models/ModelBase.h
index 4bbe31b64..d118e1516 100644
--- a/src/models/ModelBase.h
+++ b/src/models/ModelBase.h
@@ -89,6 +89,14 @@ namespace storm {
              */
             virtual bool isSymbolicModel() const;
             
+            /*!
+             * Checks whether the model is of the given type.
+             *
+             * @param modelType The model type to check for.
+             * @return True iff the model is of the given type.
+             */
+            bool isOfType(storm::models::ModelType const& modelType) const;
+            
         private:
             // The type of the model.
             ModelType modelType;
diff --git a/src/models/sparse/Model.h b/src/models/sparse/Model.h
index ce5a85b89..3d7776800 100644
--- a/src/models/sparse/Model.h
+++ b/src/models/sparse/Model.h
@@ -154,7 +154,6 @@ namespace storm {
                  */
                 RewardModelType const& getRewardModel(std::string const& rewardModelName) const;
 
-
                 /*!
                  * Retrieves the unique reward model, if there exists exactly one. Otherwise, an exception is thrown.
                  *
diff --git a/src/models/sparse/StandardRewardModel.cpp b/src/models/sparse/StandardRewardModel.cpp
index 1f7db2b79..4cab9bbb4 100644
--- a/src/models/sparse/StandardRewardModel.cpp
+++ b/src/models/sparse/StandardRewardModel.cpp
@@ -147,6 +147,7 @@ namespace storm {
                 if (this->hasTransitionRewards()) {
                     if (this->hasStateActionRewards()) {
                         storm::utility::vector::addVectors<ValueType>(this->getStateActionRewardVector(), transitionMatrix.getPointwiseProductRowSumVector(this->getTransitionRewardMatrix()), this->getStateActionRewardVector());
+                        this->optionalStateActionRewardVector = boost::none;
                     } else {
                         this->optionalStateActionRewardVector = transitionMatrix.getPointwiseProductRowSumVector(this->getTransitionRewardMatrix());
                     }
@@ -157,8 +158,9 @@ namespace storm {
                         STORM_LOG_THROW(this->getStateRewardVector().size() == this->getStateActionRewardVector().size(), storm::exceptions::InvalidOperationException, "The reduction to state rewards is only possible of both the state and the state-action rewards have the same dimension.");
                         storm::utility::vector::addVectors<ValueType>(this->getStateActionRewardVector(), this->getStateRewardVector(), this->getStateRewardVector());
                     } else {
-                        this->optionalStateRewardVector = std::move(this->optionalStateRewardVector);
+                        this->optionalStateRewardVector = std::move(this->optionalStateActionRewardVector);
                     }
+                    this->optionalStateActionRewardVector = boost::none;
                 }
             }
             
diff --git a/src/models/symbolic/Model.h b/src/models/symbolic/Model.h
index b70c39b97..4dc6ff31e 100644
--- a/src/models/symbolic/Model.h
+++ b/src/models/symbolic/Model.h
@@ -38,6 +38,7 @@ namespace storm {
             template<storm::dd::DdType Type>
             class Model : public storm::models::ModelBase {
             public:
+                static const storm::dd::DdType DdType = Type;
                 typedef StandardRewardModel<Type, double> RewardModelType;
                 
                 Model(Model<Type> const& other) = default;
diff --git a/src/settings/OptionBuilder.h b/src/settings/OptionBuilder.h
index 622a7cb1c..f70ac7c2a 100644
--- a/src/settings/OptionBuilder.h
+++ b/src/settings/OptionBuilder.h
@@ -85,7 +85,7 @@ namespace storm {
              * @return The resulting option.
              */
             std::shared_ptr<Option> build() {
-                STORM_LOG_THROW(!this->isBuild, storm::exceptions::IllegalFunctionCallException, "Cannot rebuild an option with one builder.")
+                STORM_LOG_THROW(!this->isBuild, storm::exceptions::IllegalFunctionCallException, "Cannot rebuild an option with one builder.");
 				this->isBuild = true;
 
                 if (this->hasShortName) {
diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
index 86b8dfa23..e24019e79 100644
--- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
@@ -124,7 +124,7 @@ namespace storm {
 					storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(this->policy, true);
 					submatrix.convertToEquationSystem();
 					
-					GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, this->precision, this->maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, this->relative, 50);
+					GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, this->precision, this->maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, this->relative, 50);
 					storm::utility::vector::selectVectorValues<ValueType>(subB, this->policy, rowGroupIndices, b);
 
 					// Copy X since we will overwrite it
diff --git a/src/storage/DeterministicModelBisimulationDecomposition.cpp b/src/storage/DeterministicModelBisimulationDecomposition.cpp
deleted file mode 100644
index b861f8f2c..000000000
--- a/src/storage/DeterministicModelBisimulationDecomposition.cpp
+++ /dev/null
@@ -1,1538 +0,0 @@
-#include "src/storage/DeterministicModelBisimulationDecomposition.h"
-
-#include <algorithm>
-#include <unordered_map>
-#include <chrono>
-#include <iomanip>
-#include <boost/iterator/transform_iterator.hpp>
-
-#include "src/adapters/CarlAdapter.h"
-#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
-#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
-
-#include "src/models/sparse/StandardRewardModel.h"
-
-#include "src/utility/graph.h"
-#include "src/utility/constants.h"
-#include "src/utility/ConstantsComparator.h"
-#include "src/exceptions/IllegalFunctionCallException.h"
-#include "src/exceptions/InvalidOptionException.h"
-#include "src/exceptions/InvalidArgumentException.h"
-
-#include "src/settings/SettingsManager.h"
-#include "src/settings/modules/GeneralSettings.h"
-
-namespace storm {
-    namespace storage {
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Block::Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::size_t id) : next(next), prev(prev), begin(begin), end(end), markedAsSplitter(false), markedAsPredecessorBlock(false), markedPosition(begin), absorbing(false), id(id) {
-            if (next != nullptr) {
-                next->prev = this;
-            }
-            if (prev != nullptr) {
-                prev->next = this;
-            }
-            STORM_LOG_ASSERT(begin < end, "Unable to create block of illegal size.");
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::print(Partition const& partition) const {
-            std::cout << "block " << this->getId() << " with marked position " << this->getMarkedPosition() << " and original begin " << this->getOriginalBegin() << std::endl;
-            std::cout << "begin: " << this->begin << " and end: " << this->end << " (number of states: " << this->getNumberOfStates() << ")" << std::endl;
-            std::cout << "states:" << std::endl;
-            for (storm::storage::sparse::state_type index = this->begin; index < this->end; ++index) {
-                std::cout << partition.statesAndValues[index].first << " ";
-            }
-            std::cout << std::endl << "original: " << std::endl;
-            for (storm::storage::sparse::state_type index = this->getOriginalBegin(); index < this->end; ++index) {
-                std::cout << partition.statesAndValues[index].first << " ";
-            }
-            std::cout << std::endl << "values:" << std::endl;
-            for (storm::storage::sparse::state_type index = this->begin; index < this->end; ++index) {
-                std::cout << std::setprecision(3) << partition.statesAndValues[index].second << " ";
-            }
-            if (partition.keepSilentProbabilities) {
-                std::cout << std::endl << "silent:" << std::endl;
-                for (storm::storage::sparse::state_type index = this->begin; index < this->end; ++index) {
-                    std::cout << std::setprecision(3) << partition.silentProbabilities[partition.statesAndValues[index].first] << " ";
-                }
-            }
-            std::cout << std::endl;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setBegin(storm::storage::sparse::state_type begin) {
-            this->begin = begin;
-            this->markedPosition = begin;
-            STORM_LOG_ASSERT(begin < end, "Unable to resize block to illegal size.");
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setEnd(storm::storage::sparse::state_type end) {
-            this->end = end;
-            STORM_LOG_ASSERT(begin < end, "Unable to resize block to illegal size.");
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::incrementBegin() {
-            ++this->begin;
-            STORM_LOG_ASSERT(begin <= end, "Unable to resize block to illegal size.");
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type DeterministicModelBisimulationDecomposition<ValueType>::Block::getBegin() const {
-            return this->begin;
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type DeterministicModelBisimulationDecomposition<ValueType>::Block::getOriginalBegin() const {
-            if (this->hasPreviousBlock()) {
-                return this->getPreviousBlock().getEnd();
-            } else {
-                return 0;
-            }
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type DeterministicModelBisimulationDecomposition<ValueType>::Block::getEnd() const {
-            return this->end;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setIterator(iterator it) {
-            this->selfIt = it;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block::iterator DeterministicModelBisimulationDecomposition<ValueType>::Block::getIterator() const {
-            return this->selfIt;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block::iterator DeterministicModelBisimulationDecomposition<ValueType>::Block::getNextIterator() const {
-            return this->getNextBlock().getIterator();
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block::iterator DeterministicModelBisimulationDecomposition<ValueType>::Block::getPreviousIterator() const {
-            return this->getPreviousBlock().getIterator();
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Block::getNextBlock() {
-            return *this->next;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block const& DeterministicModelBisimulationDecomposition<ValueType>::Block::getNextBlock() const {
-            return *this->next;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::hasNextBlock() const {
-            return this->next != nullptr;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Block::getPreviousBlock() {
-            return *this->prev;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block* DeterministicModelBisimulationDecomposition<ValueType>::Block::getPreviousBlockPointer() {
-            return this->prev;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block* DeterministicModelBisimulationDecomposition<ValueType>::Block::getNextBlockPointer() {
-            return this->next;
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block const& DeterministicModelBisimulationDecomposition<ValueType>::Block::getPreviousBlock() const {
-            return *this->prev;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::hasPreviousBlock() const {
-            return this->prev != nullptr;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::check() const {
-            assert(this->begin < this->end);
-            assert(this->prev == nullptr || this->prev->next == this);
-            assert(this->next == nullptr || this->next->prev == this);
-            return true;
-        }
-        
-        template<typename ValueType>
-        std::size_t DeterministicModelBisimulationDecomposition<ValueType>::Block::getNumberOfStates() const {
-            // We need to take the original begin here, because the begin is temporarily moved.
-            return (this->end - this->getOriginalBegin());
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::isMarkedAsSplitter() const {
-            return this->markedAsSplitter;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::markAsSplitter() {
-            this->markedAsSplitter = true;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::unmarkAsSplitter() {
-            this->markedAsSplitter = false;
-        }
-        
-        template<typename ValueType>
-        std::size_t DeterministicModelBisimulationDecomposition<ValueType>::Block::getId() const {
-            return this->id;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setAbsorbing(bool absorbing) {
-            this->absorbing = absorbing;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::isAbsorbing() const {
-            return this->absorbing;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setRepresentativeState(storm::storage::sparse::state_type representativeState) {
-            this->representativeState = representativeState;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::hasRepresentativeState() const {
-            return static_cast<bool>(representativeState);
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type DeterministicModelBisimulationDecomposition<ValueType>::Block::getRepresentativeState() const {
-            STORM_LOG_THROW(representativeState, storm::exceptions::InvalidOperationException, "Unable to retrieve representative state for block.");
-            return representativeState.get();
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type DeterministicModelBisimulationDecomposition<ValueType>::Block::getMarkedPosition() const {
-            return this->markedPosition;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::setMarkedPosition(storm::storage::sparse::state_type position) {
-            this->markedPosition = position;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::resetMarkedPosition() {
-            this->markedPosition = this->begin;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::incrementMarkedPosition() {
-            ++this->markedPosition;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::markAsPredecessorBlock() {
-            this->markedAsPredecessorBlock = true;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Block::unmarkAsPredecessorBlock() {
-            this->markedAsPredecessorBlock = false;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Block::isMarkedAsPredecessor() const {
-            return markedAsPredecessorBlock;
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Partition::Partition(std::size_t numberOfStates, bool keepSilentProbabilities) : numberOfBlocks(0), stateToBlockMapping(numberOfStates), statesAndValues(numberOfStates), positions(numberOfStates), keepSilentProbabilities(keepSilentProbabilities), silentProbabilities() {
-            // Create the block and give it an iterator to itself.
-            typename std::list<Block>::iterator it = blocks.emplace(this->blocks.end(), 0, numberOfStates, nullptr, nullptr, numberOfBlocks++);
-            it->setIterator(it);
-            
-            // Set up the different parts of the internal structure.
-            for (storm::storage::sparse::state_type state = 0; state < numberOfStates; ++state) {
-                statesAndValues[state] = std::make_pair(state, storm::utility::zero<ValueType>());
-                positions[state] = state;
-                stateToBlockMapping[state] = &blocks.back();
-            }
-            
-            // If we are requested to store silent probabilities, we need to prepare the vector.
-            if (this->keepSilentProbabilities) {
-                silentProbabilities = std::vector<ValueType>(numberOfStates);
-            }
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Partition::Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, boost::optional<storm::storage::sparse::state_type> representativeProb1State, bool keepSilentProbabilities) : numberOfBlocks(0), stateToBlockMapping(numberOfStates), statesAndValues(numberOfStates), positions(numberOfStates), keepSilentProbabilities(keepSilentProbabilities), silentProbabilities() {
-            storm::storage::sparse::state_type position = 0;
-            Block* firstBlock = nullptr;
-            Block* secondBlock = nullptr;
-            Block* thirdBlock = nullptr;
-            if (!prob0States.empty()) {
-                typename std::list<Block>::iterator firstIt = blocks.emplace(this->blocks.end(), 0, prob0States.getNumberOfSetBits(), nullptr, nullptr, numberOfBlocks++);
-                firstBlock = &(*firstIt);
-                firstBlock->setIterator(firstIt);
-                
-                for (auto state : prob0States) {
-                    statesAndValues[position] = std::make_pair(state, storm::utility::zero<ValueType>());
-                    positions[state] = position;
-                    stateToBlockMapping[state] = firstBlock;
-                    ++position;
-                }
-                firstBlock->setAbsorbing(true);
-            }
-            
-            if (!prob1States.empty()) {
-                typename std::list<Block>::iterator secondIt = blocks.emplace(this->blocks.end(), position, position + prob1States.getNumberOfSetBits(), firstBlock, nullptr, numberOfBlocks++);
-                secondBlock = &(*secondIt);
-                secondBlock->setIterator(secondIt);
-                
-                for (auto state : prob1States) {
-                    statesAndValues[position] = std::make_pair(state, storm::utility::zero<ValueType>());
-                    positions[state] = position;
-                    stateToBlockMapping[state] = secondBlock;
-                    ++position;
-                }
-                secondBlock->setAbsorbing(true);
-                secondBlock->setRepresentativeState(representativeProb1State.get());
-            }
-            
-            storm::storage::BitVector otherStates = ~(prob0States | prob1States);
-            if (!otherStates.empty()) {
-                typename std::list<Block>::iterator thirdIt = blocks.emplace(this->blocks.end(), position, numberOfStates, secondBlock, nullptr, numberOfBlocks++);
-                thirdBlock = &(*thirdIt);
-                thirdBlock->setIterator(thirdIt);
-                
-                for (auto state : otherStates) {
-                    statesAndValues[position] = std::make_pair(state, storm::utility::zero<ValueType>());
-                    positions[state] = position;
-                    stateToBlockMapping[state] = thirdBlock;
-                    ++position;
-                }
-            }
-            
-            // If we are requested to store silent probabilities, we need to prepare the vector.
-            if (this->keepSilentProbabilities) {
-                silentProbabilities = std::vector<ValueType>(numberOfStates);
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
-            std::swap(this->statesAndValues[this->positions[state1]], this->statesAndValues[this->positions[state2]]);
-            std::swap(this->positions[state1], this->positions[state2]);
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2) {
-            storm::storage::sparse::state_type state1 = this->statesAndValues[position1].first;
-            storm::storage::sparse::state_type state2 = this->statesAndValues[position2].first;
-            
-            std::swap(this->statesAndValues[position1], this->statesAndValues[position2]);
-            
-            this->positions[state1] = position2;
-            this->positions[state2] = position1;
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getPosition(storm::storage::sparse::state_type state) const {
-            return this->positions[state];
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position) {
-            this->positions[state] = position;
-        }
-        
-        template<typename ValueType>
-        storm::storage::sparse::state_type const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getState(storm::storage::sparse::state_type position) const {
-            return this->statesAndValues[position].first;
-        }
-        
-        template<typename ValueType>
-        ValueType const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getValue(storm::storage::sparse::state_type state) const {
-            return this->statesAndValues[this->positions[state]].second;
-        }
-        
-        template<typename ValueType>
-        ValueType const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getValueAtPosition(storm::storage::sparse::state_type position) const {
-            return this->statesAndValues[position].second;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::setValue(storm::storage::sparse::state_type state, ValueType value) {
-            this->statesAndValues[this->positions[state]].second = value;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::increaseValue(storm::storage::sparse::state_type state, ValueType value) {
-            this->statesAndValues[this->positions[state]].second += value;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::updateBlockMapping(Block& block, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last) {
-            for (; first != last; ++first) {
-                this->stateToBlockMapping[first->first] = &block;
-            }
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getFirstBlock() {
-            return *this->blocks.begin();
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBlock(storm::storage::sparse::state_type state) {
-            return *this->stateToBlockMapping[state];
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBlock(storm::storage::sparse::state_type state) const {
-            return *this->stateToBlockMapping[state];
-        }
-        
-        template<typename ValueType>
-        typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBegin(Block const& block) {
-            return this->statesAndValues.begin() + block.getBegin();
-        }
-        
-        template<typename ValueType>
-        typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator DeterministicModelBisimulationDecomposition<ValueType>::Partition::getEnd(Block const& block) {
-            return this->statesAndValues.begin() + block.getEnd();
-        }
-        
-        template<typename ValueType>
-        typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBegin(Block const& block) const {
-            return this->statesAndValues.begin() + block.getBegin();
-        }
-        
-        template<typename ValueType>
-        typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator DeterministicModelBisimulationDecomposition<ValueType>::Partition::getEnd(Block const& block) const {
-            return this->statesAndValues.begin() + block.getEnd();
-        }
-        
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Partition::splitBlock(Block& block, storm::storage::sparse::state_type position) {
-            // In case one of the resulting blocks would be empty, we simply return the current block and do not create
-            // a new one.
-            if (position == block.getBegin() || position == block.getEnd()) {
-                return block;
-            }
-            
-            // Actually create the new block and insert it at the correct position.
-            typename std::list<Block>::iterator selfIt = this->blocks.emplace(block.getIterator(), block.getBegin(), position, block.getPreviousBlockPointer(), &block, numberOfBlocks++);
-            selfIt->setIterator(selfIt);
-            Block& newBlock = *selfIt;
-            
-            // Resize the current block appropriately.
-            block.setBegin(position);
-            
-            // Update the mapping of the states in the newly created block.
-            this->updateBlockMapping(newBlock, this->getBegin(newBlock), this->getEnd(newBlock));
-            
-            return newBlock;
-        }
-
-        template<typename ValueType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Block& DeterministicModelBisimulationDecomposition<ValueType>::Partition::insertBlock(Block& block) {
-            // Find the beginning of the new block.
-            storm::storage::sparse::state_type begin;
-            if (block.hasPreviousBlock()) {
-                begin = block.getPreviousBlock().getEnd();
-            } else {
-                begin = 0;
-            }
-            
-            // Actually insert the new block.
-            typename std::list<Block>::iterator it = this->blocks.emplace(block.getIterator(), begin, block.getBegin(), block.getPreviousBlockPointer(), &block, numberOfBlocks++);
-            Block& newBlock = *it;
-            newBlock.setIterator(it);
-            
-            // Update the mapping of the states in the newly created block.
-            for (auto it = this->getBegin(newBlock), ite = this->getEnd(newBlock); it != ite; ++it) {
-                stateToBlockMapping[it->first] = &newBlock;
-            }
-            
-            return *it;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::splitLabel(storm::storage::BitVector const& statesWithLabel) {
-            for (auto blockIterator = this->blocks.begin(), ite = this->blocks.end(); blockIterator != ite; ) { // The update of the loop was intentionally moved to the bottom of the loop.
-                Block& block = *blockIterator;
-                
-                // Sort the range of the block such that all states that have the label are moved to the front.
-                std::sort(this->getBegin(block), this->getEnd(block), [&statesWithLabel] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return statesWithLabel.get(a.first) && !statesWithLabel.get(b.first); } );
-                
-                // Update the positions vector.
-                storm::storage::sparse::state_type position = block.getBegin();
-                for (auto stateIt = this->getBegin(block), stateIte = this->getEnd(block); stateIt != stateIte; ++stateIt, ++position) {
-                    this->positions[stateIt->first] = position;
-                }
-                
-                // Now we can find the first position in the block that does not have the label and create new blocks.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator it = std::find_if(this->getBegin(block), this->getEnd(block), [&] (std::pair<storm::storage::sparse::state_type, ValueType> const& a) { return !statesWithLabel.get(a.first); });
-                
-                // If not all the states agreed on the validity of the label, we need to split the block.
-                if (it != this->getBegin(block) && it != this->getEnd(block)) {
-                    auto cutPoint = std::distance(this->statesAndValues.begin(), it);
-                    this->splitBlock(block, cutPoint);
-                } else {
-                    // Otherwise, we simply proceed to the next block.
-                    ++blockIterator;
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Partition::isSilent(storm::storage::sparse::state_type state, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to retrieve silentness of state, because silent probabilities are not tracked.");
-            return comparator.isOne(this->silentProbabilities[state]);
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Partition::hasSilentProbability(storm::storage::sparse::state_type state, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to retrieve silentness of state, because silent probabilities are not tracked.");
-            return !comparator.isZero(this->silentProbabilities[state]);
-        }
-        
-        template<typename ValueType>
-        ValueType const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getSilentProbability(storm::storage::sparse::state_type state) const {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to retrieve silent probability of state, because silent probabilities are not tracked.");
-            return this->silentProbabilities[state];
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::setSilentProbabilities(typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last) {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to set silent probability of state, because silent probabilities are not tracked.");
-            for (; first != last; ++first) {
-                this->silentProbabilities[first->first] = first->second;
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::setSilentProbabilitiesToZero(typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last) {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to set silent probability of state, because silent probabilities are not tracked.");
-            for (; first != last; ++first) {
-                this->silentProbabilities[first->first] = storm::utility::zero<ValueType>();
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::setSilentProbability(storm::storage::sparse::state_type state, ValueType const& value) {
-            STORM_LOG_ASSERT(this->keepSilentProbabilities, "Unable to set silent probability of state, because silent probabilities are not tracked.");
-            this->silentProbabilities[state] = value;
-        }
-        
-        template<typename ValueType>
-        std::list<typename DeterministicModelBisimulationDecomposition<ValueType>::Block> const& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBlocks() const {
-            return this->blocks;
-        }
-        
-        template<typename ValueType>
-        std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getStatesAndValues() {
-            return this->statesAndValues;
-        }
-        
-        template<typename ValueType>
-        std::list<typename DeterministicModelBisimulationDecomposition<ValueType>::Block>& DeterministicModelBisimulationDecomposition<ValueType>::Partition::getBlocks() {
-            return this->blocks;
-        }
-        
-        template<typename ValueType>
-        bool DeterministicModelBisimulationDecomposition<ValueType>::Partition::check() const {
-            for (uint_fast64_t state = 0; state < this->positions.size(); ++state) {
-                if (this->statesAndValues[this->positions[state]].first != state) {
-                    assert(false);
-                }
-            }
-            for (auto const& block : this->blocks) {
-                assert(block.check());
-                for (auto stateIt = this->getBegin(block), stateIte = this->getEnd(block); stateIt != stateIte; ++stateIt) {
-                    if (this->stateToBlockMapping[stateIt->first] != &block) {
-                        assert(false);
-                    }
-                }
-            }
-            return true;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Partition::print() const {
-            for (auto const& block : this->blocks) {
-                block.print(*this);
-            }
-            std::cout << "states in partition" << std::endl;
-            for (auto const& stateValue : statesAndValues) {
-                std::cout << stateValue.first << " ";
-            }
-            std::cout << std::endl << "positions: " << std::endl;
-            for (auto const& index : positions) {
-                std::cout << index << " ";
-            }
-            std::cout << std::endl << "state to block mapping: " << std::endl;
-            for (auto const& block : stateToBlockMapping) {
-                std::cout << block << "[" << block->getId() <<"] ";
-            }
-            std::cout << std::endl;
-            if (this->keepSilentProbabilities) {
-                std::cout << "silent probabilities" << std::endl;
-                for (auto const& prob : silentProbabilities) {
-                    std::cout << prob << " ";
-                }
-                std::cout << std::endl;
-            }
-            std::cout << "size: " << numberOfBlocks << std::endl;
-            assert(this->check());
-        }
-        
-        template<typename ValueType>
-        std::size_t DeterministicModelBisimulationDecomposition<ValueType>::Partition::size() const {
-            return numberOfBlocks;
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Options::Options(storm::models::sparse::Model<ValueType> const& model, storm::logic::Formula const& formula) : Options() {
-            this->preserveSingleFormula(model, formula);
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Options::Options(storm::models::sparse::Model<ValueType> const& model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) : Options() {
-            if (formulas.size() == 1) {
-                this->preserveSingleFormula(model, *formulas.front());
-            } else {
-                for (auto const& formula : formulas) {
-                    std::vector<std::shared_ptr<storm::logic::AtomicExpressionFormula const>> atomicExpressionFormulas = formula->getAtomicExpressionFormulas();
-                    std::vector<std::shared_ptr<storm::logic::AtomicLabelFormula const>> atomicLabelFormulas = formula->getAtomicLabelFormulas();
-                    
-                    std::set<std::string> labelsToRespect;
-                    for (auto const& labelFormula : atomicLabelFormulas) {
-                        labelsToRespect.insert(labelFormula->getLabel());
-                    }
-                    for (auto const& expressionFormula : atomicExpressionFormulas) {
-                        std::stringstream stream;
-                        stream << *expressionFormula;
-                        labelsToRespect.insert(stream.str());
-                    }
-                    if (!respectedAtomicPropositions) {
-                        respectedAtomicPropositions = labelsToRespect;
-                    } else {
-                        respectedAtomicPropositions.get().insert(labelsToRespect.begin(), labelsToRespect.end());
-                    }
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::Options::Options() : measureDrivenInitialPartition(false), phiStates(), psiStates(), respectedAtomicPropositions(), keepRewards(true), weak(false), bounded(true), buildQuotient(true) {
-            // Intentionally left empty.
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::Options::preserveSingleFormula(storm::models::sparse::Model<ValueType> const& model, storm::logic::Formula const& formula) {
-            if (!formula.containsRewardOperator()) {
-                this->keepRewards = false;
-            }
-            
-            // As a tradeoff, we compute a strong bisimulation, because that is typically much faster. If the number of
-            // states is to be minimized, this should be set to weak by the calling site.
-            weak = false;
-            
-            // We need to preserve bounded properties iff the formula contains a bounded until or a next subformula.
-            bounded = formula.containsBoundedUntilFormula() || formula.containsNextFormula();
-            
-            std::vector<std::shared_ptr<storm::logic::AtomicExpressionFormula const>> atomicExpressionFormulas = formula.getAtomicExpressionFormulas();
-            std::vector<std::shared_ptr<storm::logic::AtomicLabelFormula const>> atomicLabelFormulas = formula.getAtomicLabelFormulas();
-            
-            std::set<std::string> labelsToRespect;
-            for (auto const& labelFormula : atomicLabelFormulas) {
-                labelsToRespect.insert(labelFormula->getLabel());
-            }
-            for (auto const& expressionFormula : atomicExpressionFormulas) {
-                std::stringstream stream;
-                stream << *expressionFormula;
-                labelsToRespect.insert(stream.str());
-            }
-            respectedAtomicPropositions = labelsToRespect;
-            
-            std::shared_ptr<storm::logic::Formula const> newFormula = formula.asSharedPointer();
-            
-            if (formula.isProbabilityOperatorFormula()) {
-                newFormula = formula.asProbabilityOperatorFormula().getSubformula().asSharedPointer();
-            } else if (formula.isRewardOperatorFormula()) {
-                newFormula = formula.asRewardOperatorFormula().getSubformula().asSharedPointer();
-            }
-            
-            std::shared_ptr<storm::logic::Formula const> leftSubformula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
-            std::shared_ptr<storm::logic::Formula const> rightSubformula;
-            if (newFormula->isUntilFormula()) {
-                leftSubformula = newFormula->asUntilFormula().getLeftSubformula().asSharedPointer();
-                rightSubformula = newFormula->asUntilFormula().getRightSubformula().asSharedPointer();
-                if (leftSubformula->isPropositionalFormula() && rightSubformula->isPropositionalFormula()) {
-                    measureDrivenInitialPartition = true;
-                }
-            } else if (newFormula->isEventuallyFormula()) {
-                rightSubformula = newFormula->asEventuallyFormula().getSubformula().asSharedPointer();
-                if (rightSubformula->isPropositionalFormula()) {
-                    measureDrivenInitialPartition = true;
-                }
-            } else if (newFormula->isReachabilityRewardFormula()) {
-                rightSubformula = newFormula->asReachabilityRewardFormula().getSubformula().asSharedPointer();
-                if (rightSubformula->isPropositionalFormula()) {
-                    measureDrivenInitialPartition = true;
-                }
-            }
-            
-            if (measureDrivenInitialPartition) {
-                storm::modelchecker::SparsePropositionalModelChecker<storm::models::sparse::Model<ValueType>> checker(model);
-                std::unique_ptr<storm::modelchecker::CheckResult> phiStatesCheckResult = checker.check(*leftSubformula);
-                std::unique_ptr<storm::modelchecker::CheckResult> psiStatesCheckResult = checker.check(*rightSubformula);
-                phiStates = phiStatesCheckResult->asExplicitQualitativeCheckResult().getTruthValuesVector();
-                psiStates = psiStatesCheckResult->asExplicitQualitativeCheckResult().getTruthValuesVector();
-            }
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::DeterministicModelBisimulationDecomposition(storm::models::sparse::Dtmc<ValueType> const& model, Options const& options) {
-            STORM_LOG_THROW(!model.hasRewardModel() || model.hasUniqueRewardModel(), storm::exceptions::IllegalFunctionCallException, "Bisimulation currently only supports models with at most one reward model.");
-            STORM_LOG_THROW(!model.hasRewardModel() || model.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently supported for models with state rewards only. Consider converting the transition rewards to state rewards (via suitable function calls).");
-            STORM_LOG_THROW(!options.weak || !options.bounded, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation cannot preserve bounded properties.");
-            storm::storage::SparseMatrix<ValueType> backwardTransitions = model.getBackwardTransitions();
-            BisimulationType bisimulationType = options.weak ? BisimulationType::WeakDtmc : BisimulationType::Strong;
-            
-            std::set<std::string> atomicPropositions;
-            if (options.respectedAtomicPropositions) {
-                atomicPropositions = options.respectedAtomicPropositions.get();
-            } else {
-                atomicPropositions = model.getStateLabeling().getLabels();
-            }
-            
-            Partition initialPartition;
-            storm::utility::ConstantsComparator<ValueType> comparator;
-            if (options.measureDrivenInitialPartition) {
-                STORM_LOG_THROW(options.phiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without phi states.");
-                STORM_LOG_THROW(options.psiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without psi states.");
-                initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, options.phiStates.get(), options.psiStates.get(), bisimulationType, options.keepRewards, options.bounded, comparator);
-            } else {
-                initialPartition = getLabelBasedInitialPartition(model, backwardTransitions, bisimulationType, atomicPropositions, options.keepRewards, comparator);
-            }
-            
-            partitionRefinement(model, atomicPropositions, backwardTransitions, initialPartition, bisimulationType, options.keepRewards, options.buildQuotient, comparator);
-        }
-        
-        template<typename ValueType>
-        DeterministicModelBisimulationDecomposition<ValueType>::DeterministicModelBisimulationDecomposition(storm::models::sparse::Ctmc<ValueType> const& model, Options const& options) {
-            STORM_LOG_THROW(!model.hasRewardModel() || model.hasUniqueRewardModel(), storm::exceptions::IllegalFunctionCallException, "Bisimulation currently only supports models with at most one reward model.");
-            STORM_LOG_THROW(!model.hasRewardModel() || model.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently supported for models with state rewards only. Consider converting the transition rewards to state rewards (via suitable function calls).");
-            STORM_LOG_THROW(!options.weak || !options.bounded, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation cannot preserve bounded properties.");
-            storm::storage::SparseMatrix<ValueType> backwardTransitions = model.getBackwardTransitions();
-            BisimulationType bisimulationType = options.weak ? BisimulationType::WeakCtmc : BisimulationType::Strong;
-            
-            std::set<std::string> atomicPropositions;
-            if (options.respectedAtomicPropositions) {
-                atomicPropositions = options.respectedAtomicPropositions.get();
-            } else {
-                atomicPropositions = model.getStateLabeling().getLabels();
-            }
-            
-            Partition initialPartition;
-            storm::utility::ConstantsComparator<ValueType> comparator;
-            if (options.measureDrivenInitialPartition) {
-                STORM_LOG_THROW(options.phiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without phi and psi states.");
-                STORM_LOG_THROW(options.psiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without phi and psi states.");
-                initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, options.phiStates.get(), options.psiStates.get(), bisimulationType, options.keepRewards, options.bounded, comparator);
-            } else {
-                initialPartition = getLabelBasedInitialPartition(model, backwardTransitions, bisimulationType, atomicPropositions, options.keepRewards, comparator);
-            }
-            
-            partitionRefinement(model, atomicPropositions, backwardTransitions, initialPartition, bisimulationType, options.keepRewards, options.buildQuotient, comparator);
-        }
-        
-        template<typename ValueType>
-        std::shared_ptr<storm::models::sparse::DeterministicModel<ValueType>> DeterministicModelBisimulationDecomposition<ValueType>::getQuotient() const {
-            STORM_LOG_THROW(this->quotient != nullptr, storm::exceptions::IllegalFunctionCallException, "Unable to retrieve quotient model from bisimulation decomposition, because it was not built.");
-            return this->quotient;
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::buildQuotient(ModelType const& model, std::set<std::string> const& selectedAtomicPropositions, Partition const& partition, BisimulationType bisimulationType, bool keepRewards, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            // In order to create the quotient model, we need to construct
-            // (a) the new transition matrix,
-            // (b) the new labeling,
-            // (c) the new reward structures.
-            
-            // Prepare a matrix builder for (a).
-            storm::storage::SparseMatrixBuilder<ValueType> builder(this->size(), this->size());
-            
-            // Prepare the new state labeling for (b).
-            storm::models::sparse::StateLabeling newLabeling(this->size());
-            std::set<std::string> atomicPropositionsSet = selectedAtomicPropositions;
-            atomicPropositionsSet.insert("init");
-            std::vector<std::string> atomicPropositions = std::vector<std::string>(atomicPropositionsSet.begin(), atomicPropositionsSet.end());
-            for (auto const& ap : atomicPropositions) {
-                newLabeling.addLabel(ap);
-            }
-            
-            // If the model had state rewards, we need to build the state rewards for the quotient as well.
-            boost::optional<std::vector<ValueType>> stateRewards;
-            if (keepRewards && model.hasRewardModel()) {
-                stateRewards = std::vector<ValueType>(this->blocks.size());
-            }
-            
-            // Now build (a) and (b) by traversing all blocks.
-            for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) {
-                auto const& block = this->blocks[blockIndex];
-                
-                // Pick one representative state. For strong bisimulation it doesn't matter which state it is, because
-                // they all behave equally.
-                storm::storage::sparse::state_type representativeState = *block.begin();
-                
-                // However, for weak bisimulation, we need to make sure the representative state is a non-silent one (if
-                // there is any such state).
-                if (bisimulationType == BisimulationType::WeakDtmc) {
-                    for (auto const& state : block) {
-                        if (!partition.isSilent(state, comparator)) {
-                            representativeState = state;
-                            break;
-                        }
-                    }
-                }
-                
-                Block const& oldBlock = partition.getBlock(representativeState);
-                
-                // If the block is absorbing, we simply add a self-loop.
-                if (oldBlock.isAbsorbing()) {
-                    builder.addNextValue(blockIndex, blockIndex, storm::utility::one<ValueType>());
-                    
-                    // If the block has a special representative state, we retrieve it now.
-                    if (oldBlock.hasRepresentativeState()) {
-                        representativeState = oldBlock.getRepresentativeState();
-                    }
-
-                    // Add all of the selected atomic propositions that hold in the representative state to the state
-                    // representing the block.
-                    for (auto const& ap : atomicPropositions) {
-                        if (model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
-                            newLabeling.addLabelToState(ap, blockIndex);
-                        }
-                    }
-                } else {
-                    // Compute the outgoing transitions of the block.
-                    std::map<storm::storage::sparse::state_type, ValueType> blockProbability;
-                    for (auto const& entry : model.getTransitionMatrix().getRow(representativeState)) {
-                        storm::storage::sparse::state_type targetBlock = partition.getBlock(entry.getColumn()).getId();
-                        
-                        // If we are computing a weak bisimulation quotient, there is no need to add self-loops.
-                        if ((bisimulationType == BisimulationType::WeakDtmc || bisimulationType == BisimulationType::WeakCtmc) && targetBlock == blockIndex) {
-                            continue;
-                        }
-                        
-                        auto probIterator = blockProbability.find(targetBlock);
-                        if (probIterator != blockProbability.end()) {
-                            probIterator->second += entry.getValue();
-                        } else {
-                            blockProbability[targetBlock] = entry.getValue();
-                        }
-                    }
-                    
-                    // Now add them to the actual matrix.
-                    for (auto const& probabilityEntry : blockProbability) {
-                        if (bisimulationType == BisimulationType::WeakDtmc) {
-                            builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second / (storm::utility::one<ValueType>() - partition.getSilentProbability(representativeState)));
-                        } else {
-                            builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second);
-                        }
-                    }
-                    
-                    // Otherwise add all atomic propositions to the equivalence class that the representative state
-                    // satisfies.
-                    for (auto const& ap : atomicPropositions) {
-                        if (model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
-                            newLabeling.addLabelToState(ap, blockIndex);
-                        }
-                    }
-                }
-                
-                // If the model has state rewards, we simply copy the state reward of the representative state, because
-                // all states in a block are guaranteed to have the same state reward.
-                if (keepRewards && model.hasRewardModel()) {
-                    typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = model.getUniqueRewardModel();
-                    stateRewards.get()[blockIndex] = nameRewardModelPair->second.getStateRewardVector()[representativeState];
-                }
-            }
-            
-            // Now check which of the blocks of the partition contain at least one initial state.
-            for (auto initialState : model.getInitialStates()) {
-                Block const& initialBlock = partition.getBlock(initialState);
-                newLabeling.addLabelToState("init", initialBlock.getId());
-            }
-            
-            // Construct the reward model mapping.
-            std::unordered_map<std::string, typename ModelType::RewardModelType> rewardModels;
-            if (keepRewards && model.hasRewardModel()) {
-                typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = model.getUniqueRewardModel();
-                rewardModels.insert(std::make_pair(nameRewardModelPair->first, typename ModelType::RewardModelType(stateRewards)));
-            }
-            
-            // Finally construct the quotient model.
-            this->quotient = std::shared_ptr<storm::models::sparse::DeterministicModel<ValueType, typename ModelType::RewardModelType>>(new ModelType(builder.build(), std::move(newLabeling), std::move(rewardModels)));
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::partitionRefinement(ModelType const& model, std::set<std::string> const& atomicPropositions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition, BisimulationType bisimulationType, bool keepRewards, bool buildQuotient, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
-            
-            // Initially, all blocks are potential splitter, so we insert them in the splitterQueue.
-            std::chrono::high_resolution_clock::time_point refinementStart = std::chrono::high_resolution_clock::now();
-            std::deque<Block*> splitterQueue;
-            std::for_each(partition.getBlocks().begin(), partition.getBlocks().end(), [&] (Block& a) { splitterQueue.push_back(&a); a.markAsSplitter(); });
-            
-            // Then perform the actual splitting until there are no more splitters.
-            while (!splitterQueue.empty()) {
-                // Optionally: sort the splitter queue according to some criterion (here: prefer small splitters).
-                std::sort(splitterQueue.begin(), splitterQueue.end(), [] (Block const* b1, Block const* b2) { return b1->getNumberOfStates() < b2->getNumberOfStates() || (b1->getNumberOfStates() == b2->getNumberOfStates() && b1->getId() < b2->getId()); } );
-                std::sort(splitterQueue.begin(), splitterQueue.end(), [] (Block const* b1, Block const* b2) { return b1->getNumberOfStates() < b2->getNumberOfStates(); } );
-
-                // Get and prepare the next splitter.
-                Block* splitter = splitterQueue.front();
-                splitterQueue.pop_front();
-                splitter->unmarkAsSplitter();
-                
-                // Now refine the partition using the current splitter.
-                refinePartition(model.getTransitionMatrix(), backwardTransitions, *splitter, partition, bisimulationType, splitterQueue, comparator);
-            }
-            std::chrono::high_resolution_clock::duration refinementTime = std::chrono::high_resolution_clock::now() - refinementStart;
-            
-            // Now move the states from the internal partition into their final place in the decomposition. We do so in
-            // a way that maintains the block IDs as indices.
-            std::chrono::high_resolution_clock::time_point extractionStart = std::chrono::high_resolution_clock::now();
-            this->blocks.resize(partition.size());
-            for (auto const& block : partition.getBlocks()) {
-                // We need to sort the states to allow for rapid construction of the blocks.
-                std::sort(partition.getBegin(block), partition.getEnd(block), [] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return a.first < b.first; });
-                
-                // Convert the state-value-pairs to states only.
-				std::function<storm::storage::sparse::state_type (std::pair<storm::storage::sparse::state_type, ValueType> const&)> projection = [] (std::pair<storm::storage::sparse::state_type, ValueType> const& a) -> storm::storage::sparse::state_type { return a.first; };
-				this->blocks[block.getId()] = block_type(boost::make_transform_iterator(partition.getBegin(block), projection), boost::make_transform_iterator(partition.getEnd(block), projection), true);
-            }
-            
-            // If we are required to build the quotient model, do so now.
-            if (buildQuotient) {
-                this->buildQuotient(model, atomicPropositions, partition, bisimulationType, keepRewards, comparator);
-            }
-            
-            std::chrono::high_resolution_clock::duration extractionTime = std::chrono::high_resolution_clock::now() - extractionStart;
-            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
-            
-            if (storm::settings::generalSettings().isShowStatisticsSet()) {
-                std::chrono::milliseconds refinementTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(refinementTime);
-                std::chrono::milliseconds extractionTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(extractionTime);
-                std::chrono::milliseconds totalTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(totalTime);
-                std::cout << std::endl;
-                std::cout << "Time breakdown:" << std::endl;
-                std::cout << "    * time for partitioning: " << refinementTimeInMilliseconds.count() << "ms" << std::endl;
-                std::cout << "    * time for extraction: " << extractionTimeInMilliseconds.count() << "ms" << std::endl;
-                std::cout << "------------------------------------------" << std::endl;
-                std::cout << "    * total time: " << totalTimeInMilliseconds.count() << "ms" << std::endl;
-                std::cout << std::endl;
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::refineBlockProbabilities(Block& block, Partition& partition, BisimulationType bisimulationType, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            // Sort the states in the block based on their probabilities.
-            std::sort(partition.getBegin(block), partition.getEnd(block), [&partition] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return a.second < b.second; } );
-            
-            // Update the positions vector.
-            storm::storage::sparse::state_type position = block.getBegin();
-            for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt, ++position) {
-                partition.setPosition(stateIt->first, position);
-            }
-            
-            // Finally, we need to scan the ranges of states that agree on the probability.
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator begin = partition.getBegin(block);
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator end = partition.getEnd(block) - 1;
-            storm::storage::sparse::state_type currentIndex = block.getBegin();
-            
-            // Now we can check whether the block needs to be split, which is the case iff the probabilities for the
-            // first and the last state are different.
-            bool blockSplit = !comparator.isEqual(begin->second, end->second);
-            while (!comparator.isEqual(begin->second, end->second)) {
-                // Now we scan for the first state in the block that disagrees on the probability value.
-                // Note that we do not have to check currentIndex for staying within bounds, because we know the matching
-                // state is within bounds.
-                ValueType const& currentValue = begin->second;
-                
-                ++begin;
-                ++currentIndex;
-                while (begin != end && comparator.isEqual(begin->second, currentValue)) {
-                    ++begin;
-                    ++currentIndex;
-                }
-                
-                // Now we split the block and mark it as a potential splitter.
-                Block& newBlock = partition.splitBlock(block, currentIndex);
-                if (!newBlock.isMarkedAsSplitter()) {
-                    splitterQueue.push_back(&newBlock);
-                    newBlock.markAsSplitter();
-                }
-            }
-            
-            // If the block was split, we also need to insert itself into the splitter queue.
-            if (blockSplit) {
-                if (!block.isMarkedAsSplitter()) {
-                    splitterQueue.push_back(&block);
-                    block.markAsSplitter();
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::refineBlockWeak(Block& block, Partition& partition, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            std::vector<uint_fast64_t> splitPoints = getSplitPointsWeak(block, partition, comparator);
-            
-            // Restore the original begin of the block.
-            block.setBegin(block.getOriginalBegin());
-            
-            // Now that we have the split points of the non-silent states, we perform a backward search from
-            // each non-silent state and label the predecessors with the class of the non-silent state.
-            std::vector<storm::storage::BitVector> stateLabels(block.getEnd() - block.getBegin(), storm::storage::BitVector(splitPoints.size() - 1));
-            
-            std::vector<storm::storage::sparse::state_type> stateStack;
-            stateStack.reserve(block.getEnd() - block.getBegin());
-            for (uint_fast64_t stateClassIndex = 0; stateClassIndex < splitPoints.size() - 1; ++stateClassIndex) {
-                for (auto stateIt = partition.getStatesAndValues().begin() + splitPoints[stateClassIndex], stateIte = partition.getStatesAndValues().begin() + splitPoints[stateClassIndex + 1]; stateIt != stateIte; ++stateIt) {
-                    
-                    stateStack.push_back(stateIt->first);
-                    stateLabels[partition.getPosition(stateIt->first) - block.getBegin()].set(stateClassIndex);
-                    while (!stateStack.empty()) {
-                        storm::storage::sparse::state_type currentState = stateStack.back();
-                        stateStack.pop_back();
-                        
-                        for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) {
-                            if (comparator.isZero(predecessorEntry.getValue())) {
-                                continue;
-                            }
-                            
-                            storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
-                            
-                            // Only if the state is in the same block, is a silent state and it has not yet been
-                            // labeled with the current label.
-                            if (&partition.getBlock(predecessor) == &block && partition.isSilent(predecessor, comparator) && !stateLabels[partition.getPosition(predecessor) - block.getBegin()].get(stateClassIndex)) {
-                                stateStack.push_back(predecessor);
-                                stateLabels[partition.getPosition(predecessor) - block.getBegin()].set(stateClassIndex);
-                            }
-                        }
-                    }
-                }
-            }
-            
-            // Now that all states were appropriately labeled, we can sort the states according to their labels and then
-            // scan for ranges that agree on the label.
-            std::sort(partition.getBegin(block), partition.getEnd(block), [&] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return stateLabels[partition.getPosition(a.first) - block.getBegin()] < stateLabels[partition.getPosition(b.first) - block.getBegin()]; });
-            
-            // Note that we do not yet repair the positions vector, but for the sake of efficiency temporariliy keep the
-            // data structure in an inconsistent state.
-            
-            // Now we have everything in place to actually split the block by just scanning for ranges of equal label.
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator begin = partition.getBegin(block);
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator end = partition.getEnd(block) - 1;
-            storm::storage::sparse::state_type currentIndex = block.getBegin();
-            
-            // Now we can check whether the block needs to be split, which is the case iff the labels for the first and
-            // the last state are different. Store the offset of the block seperately, because it will potentially
-            // modified by splits.
-            storm::storage::sparse::state_type blockOffset = block.getBegin();
-            bool blockSplit = stateLabels[partition.getPosition(begin->first) - blockOffset] != stateLabels[partition.getPosition(end->first) - blockOffset];
-            while (stateLabels[partition.getPosition(begin->first) - blockOffset] != stateLabels[partition.getPosition(end->first) - blockOffset]) {
-                // Now we scan for the first state in the block that disagrees on the labeling value.
-                // Note that we do not have to check currentIndex for staying within bounds, because we know the matching
-                // state is within bounds.
-                storm::storage::BitVector const& currentValue = stateLabels[partition.getPosition(begin->first) - blockOffset];
-                
-                ++begin;
-                ++currentIndex;
-                while (begin != end && stateLabels[partition.getPosition(begin->first) - blockOffset] == currentValue) {
-                    ++begin;
-                    ++currentIndex;
-                }
-                
-                // Now we split the block and mark it as a potential splitter.
-                Block& newBlock = partition.splitBlock(block, currentIndex);
-                
-                // Update the silent probabilities for all the states in the new block.
-                for (auto stateIt = partition.getBegin(newBlock), stateIte = partition.getEnd(newBlock); stateIt != stateIte; ++stateIt) {
-                    if (partition.hasSilentProbability(stateIt->first, comparator)) {
-                        ValueType newSilentProbability = storm::utility::zero<ValueType>();
-                        for (auto const& successorEntry : forwardTransitions.getRow(stateIt->first)) {
-                            if (&partition.getBlock(successorEntry.getColumn()) == &newBlock) {
-                                newSilentProbability += successorEntry.getValue();
-                            }
-                        }
-                        partition.setSilentProbability(stateIt->first, newSilentProbability);
-                    }
-                }
-                
-                if (!newBlock.isMarkedAsSplitter()) {
-                    splitterQueue.push_back(&newBlock);
-                    newBlock.markAsSplitter();
-                }
-            }
-            
-            // If the block was split, we also need to insert itself into the splitter queue.
-            if (blockSplit) {
-                if (!block.isMarkedAsSplitter()) {
-                    splitterQueue.push_back(&block);
-                    block.markAsSplitter();
-                }
-                
-                // Update the silent probabilities for all the states in the old block.
-                for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt) {
-                    if (partition.hasSilentProbability(stateIt->first, comparator)) {
-                        ValueType newSilentProbability = storm::utility::zero<ValueType>();
-                        for (auto const& successorEntry : forwardTransitions.getRow(stateIt->first)) {
-                            if (&partition.getBlock(successorEntry.getColumn()) == &block) {
-                                newSilentProbability += successorEntry.getValue();
-                            }
-                        }
-                        partition.setSilentProbability(stateIt->first, newSilentProbability);
-                    }
-                }
-            }
-            
-            // Finally update the positions vector.
-            storm::storage::sparse::state_type position = blockOffset;
-            for (auto stateIt = partition.getStatesAndValues().begin() + blockOffset, stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt, ++position) {
-                partition.setPosition(stateIt->first, position);
-            }
-        }
-        
-        template<typename ValueType>
-        std::vector<uint_fast64_t> DeterministicModelBisimulationDecomposition<ValueType>::getSplitPointsWeak(Block& block, Partition& partition, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            std::vector<uint_fast64_t> result;
-            // We first scale all probabilities with (1-p[s]) where p[s] is the silent probability of state s.
-            std::for_each(partition.getStatesAndValues().begin() + block.getOriginalBegin(), partition.getStatesAndValues().begin() + block.getBegin(), [&] (std::pair<storm::storage::sparse::state_type, ValueType>& stateValuePair) {
-                ValueType const& silentProbability = partition.getSilentProbability(stateValuePair.first);
-                if (!comparator.isOne(silentProbability) && !comparator.isZero(silentProbability)) {
-                    stateValuePair.second /= storm::utility::one<ValueType>() - silentProbability;
-                }
-            });
-            
-            // Now sort the states based on their probabilities.
-            std::sort(partition.getStatesAndValues().begin() + block.getOriginalBegin(), partition.getStatesAndValues().begin() + block.getBegin(), [&partition] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return a.second < b.second; } );
-            
-            // Update the positions vector.
-            storm::storage::sparse::state_type position = block.getOriginalBegin();
-            for (auto stateIt = partition.getStatesAndValues().begin() + block.getOriginalBegin(), stateIte = partition.getStatesAndValues().begin() + block.getBegin(); stateIt != stateIte; ++stateIt, ++position) {
-                partition.setPosition(stateIt->first, position);
-            }
-            
-            // Then, we scan for the ranges of states that agree on the probability.
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator begin = partition.getStatesAndValues().begin() + block.getOriginalBegin();
-            typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator end = partition.getStatesAndValues().begin() + block.getBegin() - 1;
-            storm::storage::sparse::state_type currentIndex = block.getOriginalBegin();
-            result.push_back(currentIndex);
-            
-            // Now we can check whether the block needs to be split, which is the case iff the probabilities for the
-            // first and the last state are different.
-            while (!comparator.isEqual(begin->second, end->second)) {
-                // Now we scan for the first state in the block that disagrees on the probability value.
-                // Note that we do not have to check currentIndex for staying within bounds, because we know the matching
-                // state is within bounds.
-                ValueType const& currentValue = begin->second;
-                
-                ++begin;
-                ++currentIndex;
-                while (begin != end && comparator.isEqual(begin->second, currentValue)) {
-                    ++begin;
-                    ++currentIndex;
-                }
-                
-                // Remember the index at which the probabilities were different.
-                result.push_back(currentIndex);
-            }
-            
-            // Push a sentinel element and return result.
-            result.push_back(block.getBegin());
-            return result;
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::refinePartition(storm::storage::SparseMatrix<ValueType> const& forwardTransitions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Block& splitter, Partition& partition, BisimulationType bisimulationType, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            std::list<Block*> predecessorBlocks;
-            
-            // Iterate over all states of the splitter and check its predecessors.
-            bool splitterIsPredecessor = false;
-            storm::storage::sparse::state_type currentPosition = splitter.getBegin();
-            for (auto stateIterator = partition.getBegin(splitter), stateIte = partition.getEnd(splitter); stateIterator != stateIte; ++stateIterator, ++currentPosition) {
-                storm::storage::sparse::state_type currentState = stateIterator->first;
-                
-                uint_fast64_t elementsToSkip = 0;
-                for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) {
-                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
-                    
-                    // Get predecessor block and remember if the splitter was a predecessor of itself.
-                    Block& predecessorBlock = partition.getBlock(predecessor);
-                    if (&predecessorBlock == &splitter) {
-                        splitterIsPredecessor = true;
-                    }
-                    
-                    // If the predecessor block has just one state or is marked as being absorbing, we must not split it.
-                    if (predecessorBlock.getNumberOfStates() <= 1 || predecessorBlock.isAbsorbing()) {
-                        continue;
-                    }
-                    
-                    storm::storage::sparse::state_type predecessorPosition = partition.getPosition(predecessor);
-                    
-                    // If we have not seen this predecessor before, we move it to a part near the beginning of the block.
-                    if (predecessorPosition >= predecessorBlock.getBegin()) {
-                        if (&predecessorBlock == &splitter) {
-                            // If the predecessor we just found was already processed (in terms of visiting its predecessors),
-                            // we swap it with the state that is currently at the beginning of the block and move the start
-                            // of the block one step further.
-                            if (predecessorPosition <= currentPosition + elementsToSkip) {
-                                partition.swapStates(predecessor, partition.getState(predecessorBlock.getBegin()));
-                                predecessorBlock.incrementBegin();
-                            } else {
-                                // Otherwise, we need to move the predecessor, but we need to make sure that we explore its
-                                // predecessors later.
-                                if (predecessorBlock.getMarkedPosition() == predecessorBlock.getBegin()) {
-                                    partition.swapStatesAtPositions(predecessorBlock.getMarkedPosition(), predecessorPosition);
-                                    partition.swapStatesAtPositions(predecessorPosition, currentPosition + elementsToSkip + 1);
-                                } else {
-                                    partition.swapStatesAtPositions(predecessorBlock.getMarkedPosition(), predecessorPosition);
-                                    partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.getBegin());
-                                    partition.swapStatesAtPositions(predecessorPosition, currentPosition + elementsToSkip + 1);
-                                }
-                                
-                                ++elementsToSkip;
-                                predecessorBlock.incrementMarkedPosition();
-                                predecessorBlock.incrementBegin();
-                            }
-                        } else {
-                            partition.swapStates(predecessor, partition.getState(predecessorBlock.getBegin()));
-                            predecessorBlock.incrementBegin();
-                        }
-                        partition.setValue(predecessor, predecessorEntry.getValue());
-                    } else {
-                        // Otherwise, we just need to update the probability for this predecessor.
-                        partition.increaseValue(predecessor, predecessorEntry.getValue());
-                    }
-                    
-                    if (!predecessorBlock.isMarkedAsPredecessor()) {
-                        predecessorBlocks.emplace_back(&predecessorBlock);
-                        predecessorBlock.markAsPredecessorBlock();
-                    }
-                }
-                
-                // If we had to move some elements beyond the current element, we may have to skip them.
-                if (elementsToSkip > 0) {
-                    stateIterator += elementsToSkip;
-                    currentPosition += elementsToSkip;
-                }
-            }
-            
-            // Now we can traverse the list of states of the splitter whose predecessors we have not yet explored.
-            for (auto stateIterator = partition.getStatesAndValues().begin() + splitter.getOriginalBegin(), stateIte = partition.getStatesAndValues().begin() + splitter.getMarkedPosition(); stateIterator != stateIte; ++stateIterator) {
-                storm::storage::sparse::state_type currentState = stateIterator->first;
-                
-                for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) {
-                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
-                    Block& predecessorBlock = partition.getBlock(predecessor);
-                    storm::storage::sparse::state_type predecessorPosition = partition.getPosition(predecessor);
-                    
-                    if (predecessorPosition >= predecessorBlock.getBegin()) {
-                        partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.getBegin());
-                        predecessorBlock.incrementBegin();
-                        partition.setValue(predecessor, predecessorEntry.getValue());
-                    } else {
-                        partition.increaseValue(predecessor, predecessorEntry.getValue());
-                    }
-                    
-                    if (!predecessorBlock.isMarkedAsPredecessor()) {
-                        predecessorBlocks.emplace_back(&predecessorBlock);
-                        predecessorBlock.markAsPredecessorBlock();
-                    }
-                }
-            }
-            
-            if (bisimulationType == BisimulationType::Strong || bisimulationType == BisimulationType::WeakCtmc) {
-                std::vector<Block*> blocksToSplit;
-                
-                // Now, we can iterate over the predecessor blocks and see whether we have to create a new block for
-                // predecessors of the splitter.
-                for (auto blockPtr : predecessorBlocks) {
-                    Block& block = *blockPtr;
-                    
-                    block.unmarkAsPredecessorBlock();
-                    block.resetMarkedPosition();
-                    
-                    // If we have moved the begin of the block to somewhere in the middle of the block, we need to split it.
-                    if (block.getBegin() != block.getEnd()) {
-                        Block& newBlock = partition.insertBlock(block);
-                        if (!newBlock.isMarkedAsSplitter()) {
-                            splitterQueue.push_back(&newBlock);
-                            newBlock.markAsSplitter();
-                        }
-                        
-                        // Schedule the block of predecessors for refinement based on probabilities.
-                        blocksToSplit.emplace_back(&newBlock);
-                    } else {
-                        // In this case, we can keep the block by setting its begin to the old value.
-                        block.setBegin(block.getOriginalBegin());
-                        blocksToSplit.emplace_back(&block);
-                    }
-                }
-                
-                // Finally, we walk through the blocks that have a transition to the splitter and split them using
-                // probabilistic information.
-                for (auto blockPtr : blocksToSplit) {
-                    if (blockPtr->getNumberOfStates() <= 1) {
-                        continue;
-                    }
-                    
-                    // In the case of weak bisimulation for CTMCs, we don't need to make sure the rate of staying inside
-                    // the own block is the same.
-                    if (bisimulationType == BisimulationType::WeakCtmc && blockPtr == &splitter) {
-                        continue;
-                    }
-                    
-                    refineBlockProbabilities(*blockPtr, partition, bisimulationType, splitterQueue, comparator);
-                }
-            } else { // In this case, we are computing a weak bisimulation on a DTMC.
-                // If the splitter was a predecessor of itself and we are computing a weak bisimulation, we need to update
-                // the silent probabilities.
-                if (splitterIsPredecessor) {
-                    partition.setSilentProbabilities(partition.getStatesAndValues().begin() + splitter.getOriginalBegin(), partition.getStatesAndValues().begin() + splitter.getBegin());
-                    partition.setSilentProbabilitiesToZero(partition.getStatesAndValues().begin() + splitter.getBegin(), partition.getStatesAndValues().begin() + splitter.getEnd());
-                }
-                
-                // Now refine all predecessor blocks in a weak manner. That is, we split according to the criterion of
-                // weak bisimulation.
-                for (auto blockPtr : predecessorBlocks) {
-                    Block& block = *blockPtr;
-                    
-                    // If the splitter is also the predecessor block, we must not refine it at this point.
-                    if (&block != &splitter) {
-                        refineBlockWeak(block, partition, forwardTransitions, backwardTransitions, splitterQueue, comparator);
-                    } else {
-                        // Restore the begin of the block.
-                        block.setBegin(block.getOriginalBegin());
-                    }
-                    
-                    block.unmarkAsPredecessorBlock();
-                    block.resetMarkedPosition();
-                }
-            }
-            
-            STORM_LOG_ASSERT(partition.check(), "Partition became inconsistent.");
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Partition DeterministicModelBisimulationDecomposition<ValueType>::getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, BisimulationType bisimulationType, bool keepRewards, bool bounded, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates);
-            
-            boost::optional<storm::storage::sparse::state_type> representativePsiState;
-            if (!psiStates.empty()) {
-                representativePsiState = *psiStates.begin();
-            }
-            
-            Partition partition(model.getNumberOfStates(), statesWithProbability01.first, bounded || keepRewards ? psiStates : statesWithProbability01.second, representativePsiState, bisimulationType == BisimulationType::WeakDtmc);
-            
-            // If the model has state rewards, we need to consider them, because otherwise reward properties are not
-            // preserved.
-            if (keepRewards && model.hasRewardModel()) {
-                this->splitRewards(model.getUniqueRewardModel()->second.getStateRewardVector(), partition, comparator);
-            }
-            
-            // If we are creating the initial partition for weak bisimulation, we need to (a) split off all divergent
-            // states of each initial block and (b) initialize the vector of silent probabilities held by the partition.
-            if (bisimulationType == BisimulationType::WeakDtmc) {
-                this->splitOffDivergentStates(model, backwardTransitions, partition);
-                this->initializeSilentProbabilities(model, partition);
-            }
-            
-            return partition;
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::splitOffDivergentStates(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition) {
-            std::vector<storm::storage::sparse::state_type> stateStack;
-            stateStack.reserve(model.getNumberOfStates());
-            storm::storage::BitVector nondivergentStates(model.getNumberOfStates());
-            for (auto& block : partition.getBlocks()) {
-                nondivergentStates.clear();
-                
-                for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt) {
-                    if (nondivergentStates.get(stateIt->first)) {
-                        continue;
-                    }
-                    
-                    // Now traverse the forward transitions of the current state and check whether there is a
-                    // transition to some other block.
-                    bool isDirectlyNonDivergent = false;
-                    for (auto const& successor : model.getRows(stateIt->first)) {
-                        // If there is such a transition, then we can mark all states in the current block that can
-                        // reach the state as non-divergent.
-                        if (&partition.getBlock(successor.getColumn()) != &block) {
-                            isDirectlyNonDivergent = true;
-                            break;
-                        }
-                    }
-                    
-                    if (isDirectlyNonDivergent) {
-                        stateStack.push_back(stateIt->first);
-                        
-                        while (!stateStack.empty()) {
-                            storm::storage::sparse::state_type currentState = stateStack.back();
-                            stateStack.pop_back();
-                            nondivergentStates.set(currentState);
-                            
-                            for (auto const& predecessor : backwardTransitions.getRow(currentState)) {
-                                if (&partition.getBlock(predecessor.getColumn()) == &block && !nondivergentStates.get(predecessor.getColumn())) {
-                                    stateStack.push_back(predecessor.getColumn());
-                                }
-                            }
-                        }
-                    }
-                }
-                
-                
-                if (nondivergentStates.getNumberOfSetBits() > 0 && nondivergentStates.getNumberOfSetBits() != block.getNumberOfStates()) {
-                    // Now that we have determined all (non)divergent states in the current block, we need to split them
-                    // off.
-                    std::sort(partition.getBegin(block), partition.getEnd(block), [&nondivergentStates] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return nondivergentStates.get(a.first) && !nondivergentStates.get(b.first); } );
-                    // Update the positions vector.
-                    storm::storage::sparse::state_type position = block.getBegin();
-                    for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt, ++position) {
-                        partition.setPosition(stateIt->first, position);
-                    }
-                    
-                    // Finally, split the block.
-                    partition.splitBlock(block, block.getBegin() + nondivergentStates.getNumberOfSetBits());
-                    
-                    // Since the remaining states in the block are divergent, we can mark the block as absorbing.
-                    // This also guarantees that the self-loop will be added to the state of the quotient
-                    // representing this block of states.
-                    block.setAbsorbing(true);
-                } else if (nondivergentStates.getNumberOfSetBits() == 0) {
-                    // If there are only diverging states in the block, we need to make it absorbing.
-                    block.setAbsorbing(true);
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        typename DeterministicModelBisimulationDecomposition<ValueType>::Partition DeterministicModelBisimulationDecomposition<ValueType>::getLabelBasedInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, BisimulationType bisimulationType, boost::optional<std::set<std::string>> const& atomicPropositions, bool keepRewards, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            Partition partition(model.getNumberOfStates(), bisimulationType == BisimulationType::WeakDtmc);
-            
-            if (atomicPropositions) {
-                for (auto const& label : atomicPropositions.get()) {
-                    if (label == "init") {
-                        continue;
-                    }
-                    partition.splitLabel(model.getStates(label));
-                }
-            } else {
-                for (auto const& label : model.getStateLabeling().getLabels()) {
-                    if (label == "init") {
-                        continue;
-                    }
-                    partition.splitLabel(model.getStates(label));
-                }
-            }
-            
-            // If the model has state rewards, we need to consider them, because otherwise reward properties are not
-            // preserved.
-            if (keepRewards && model.hasRewardModel()) {
-                this->splitRewards(model.getUniqueRewardModel()->second.getStateRewardVector(), partition, comparator);
-            }
-            
-            // If we are creating the initial partition for weak bisimulation, we need to (a) split off all divergent
-            // states of each initial block and (b) initialize the vector of silent probabilities held by the partition.
-            if (bisimulationType == BisimulationType::WeakDtmc) {
-                this->splitOffDivergentStates(model, backwardTransitions, partition);
-                this->initializeSilentProbabilities(model, partition);
-            }
-            return partition;
-        }
-        
-        template<typename ValueType>
-        template<typename ModelType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::initializeSilentProbabilities(ModelType const& model, Partition& partition) {
-            for (auto const& block : partition.getBlocks()) {
-                for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt) {
-                    ValueType silentProbability = storm::utility::zero<ValueType>();
-                    
-                    for (auto const& successorEntry : model.getRows(stateIt->first)) {
-                        if (&partition.getBlock(successorEntry.getColumn()) == &block) {
-                            silentProbability += successorEntry.getValue();
-                        }
-                    }
-                    
-                    partition.setSilentProbability(stateIt->first, silentProbability);
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        void DeterministicModelBisimulationDecomposition<ValueType>::splitRewards(std::vector<ValueType> const& stateRewardVector, Partition& partition, storm::utility::ConstantsComparator<ValueType> const& comparator) {
-            for (auto& block : partition.getBlocks()) {
-                std::sort(partition.getBegin(block), partition.getEnd(block), [&stateRewardVector] (std::pair<storm::storage::sparse::state_type, ValueType> const& a, std::pair<storm::storage::sparse::state_type, ValueType> const& b) { return stateRewardVector[a.first] < stateRewardVector[b.first]; } );
-                
-                // Update the positions vector and put the (state) reward values next to the states so we can easily compare them later.
-                storm::storage::sparse::state_type position = block.getBegin();
-                for (auto stateIt = partition.getBegin(block), stateIte = partition.getEnd(block); stateIt != stateIte; ++stateIt, ++position) {
-                    partition.setPosition(stateIt->first, position);
-                    stateIt->second = stateRewardVector[stateIt->first];
-                }
-                
-                // Finally, we need to scan the ranges of states that agree on the probability.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator begin = partition.getBegin(block);
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator end = partition.getEnd(block) - 1;
-                storm::storage::sparse::state_type currentIndex = block.getBegin();
-                
-                // Now we can check whether the block needs to be split, which is the case iff the rewards for the first
-                // and the last state are different.
-                while (!comparator.isEqual(begin->second, end->second)) {
-                    // Now we scan for the first state in the block that disagrees on the reward value. Note that we do
-                    // not have to check currentIndex for staying within bounds, because we know the matching state is
-                    // within bounds.
-                    ValueType const& currentValue = begin->second;
-                    
-                    ++begin;
-                    ++currentIndex;
-                    while (begin != end && comparator.isEqual(begin->second, currentValue)) {
-                        ++begin;
-                        ++currentIndex;
-                    }
-                    
-                    // Now we split the block.
-                    partition.splitBlock(block, currentIndex);
-                }
-            }
-        }
-        
-        template class DeterministicModelBisimulationDecomposition<double>;
-        
-#ifdef STORM_HAVE_CARL
-        template class DeterministicModelBisimulationDecomposition<storm::RationalFunction>;
-#endif
-    }
-}
diff --git a/src/storage/DeterministicModelBisimulationDecomposition.h b/src/storage/DeterministicModelBisimulationDecomposition.h
deleted file mode 100644
index d66591f76..000000000
--- a/src/storage/DeterministicModelBisimulationDecomposition.h
+++ /dev/null
@@ -1,590 +0,0 @@
-#ifndef STORM_STORAGE_DeterministicModelBisimulationDecomposition_H_
-#define STORM_STORAGE_DeterministicModelBisimulationDecomposition_H_
-
-#include <queue>
-#include <deque>
-
-#include "src/storage/sparse/StateType.h"
-#include "src/storage/Decomposition.h"
-#include "src/storage/StateBlock.h"
-#include "src/logic/Formulas.h"
-#include "src/models/sparse/Dtmc.h"
-#include "src/models/sparse/Ctmc.h"
-#include "src/storage/Distribution.h"
-#include "src/utility/constants.h"
-#include "src/utility/OsDetection.h"
-#include "src/exceptions/InvalidOperationException.h"
-
-namespace storm {
-    namespace utility {
-        template <typename ValueType> class ConstantsComparator;
-    }
-    
-    namespace storage {
-        
-        /*!
-         * This class represents the decomposition model into its (strong) bisimulation quotient.
-         */
-        template <typename ValueType>
-        class DeterministicModelBisimulationDecomposition : public Decomposition<StateBlock> {
-        public:
-            // A class that offers the possibility to customize the bisimulation.
-            struct Options {
-                // Creates an object representing the default values for all options.
-                Options();
-
-                /*!
-                 * Creates an object representing the options necessary to obtain the quotient while still preserving
-                 * the given formula.
-                 *
-                 * @param The model for which the quotient model shall be computed. This needs to be given in order to
-                 * derive a suitable initial partition.
-                 * @param formula The formula that is to be preserved.
-                 */
-                Options(storm::models::sparse::Model<ValueType> const& model, storm::logic::Formula const& formula);
-
-                /*!
-                 * Creates an object representing the options necessary to obtain the smallest quotient while still
-                 * preserving the given formula.
-                 *
-                 * @param The model for which the quotient model shall be computed. This needs to be given in order to
-                 * derive a suitable initial partition.
-                 * @param formulas The formula that is to be preserved.
-                 */
-                Options(storm::models::sparse::Model<ValueType> const& model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas);
-                
-                // A flag that indicates whether a measure driven initial partition is to be used. If this flag is set
-                // to true, the two optional pairs phiStatesAndLabel and psiStatesAndLabel must be set. Then, the
-                // measure driven initial partition wrt. to the states phi and psi is taken.
-                bool measureDrivenInitialPartition;
-                boost::optional<storm::storage::BitVector> phiStates;
-                boost::optional<storm::storage::BitVector> psiStates;
-                
-                // An optional set of strings that indicate which of the atomic propositions of the model are to be
-                // respected and which may be ignored. If not given, all atomic propositions of the model are respected.
-                boost::optional<std::set<std::string>> respectedAtomicPropositions;
-                
-                // A flag that indicates whether or not the state-rewards of the model are to be respected (and should
-                // be kept in the quotient model, if one is built).
-                bool keepRewards;
-                
-                // A flag that indicates whether a strong or a weak bisimulation is to be computed.
-                bool weak;
-                
-                // A flag that indicates whether step-bounded properties are to be preserved. This may only be set to tru
-                // when computing strong bisimulation equivalence.
-                bool bounded;
-                
-                // A flag that governs whether the quotient model is actually built or only the decomposition is computed.
-                bool buildQuotient;
-                
-            private:
-                /*!
-                 * Sets the options under the assumption that the given formula is the only one that is to be checked.
-                 *
-                 * @param model The model for which to preserve the formula.
-                 * @param formula The only formula to check.
-                 */
-                void preserveSingleFormula(storm::models::sparse::Model<ValueType> const& model, storm::logic::Formula const& formula);
-            };
-            
-            /*!
-             * Decomposes the given DTMC into equivalance classes of a bisimulation. Which kind of bisimulation is
-             * computed, is customizable via the given options.
-             *
-             * @param model The model to decompose.
-             * @param options The options that customize the computed bisimulation.
-             */
-            DeterministicModelBisimulationDecomposition(storm::models::sparse::Dtmc<ValueType> const& model, Options const& options = Options());
-            
-            /*!
-             * Decomposes the given CTMC into equivalance classes of a bisimulation. Which kind of bisimulation is
-             * computed, is customizable via the given options.
-             *
-             * @param model The model to decompose.
-             * @param options The options that customize the computed bisimulation.
-             */
-            DeterministicModelBisimulationDecomposition(storm::models::sparse::Ctmc<ValueType> const& model, Options const& options = Options());
-            
-            /*!
-             * Retrieves the quotient of the model under the previously computed bisimulation.
-             *
-             * @return The quotient model.
-             */
-            std::shared_ptr<storm::models::sparse::DeterministicModel<ValueType>> getQuotient() const;
-            
-        private:
-            enum class BisimulationType { Strong, WeakDtmc, WeakCtmc };
-            
-            class Partition;
-            
-            class Block {
-            public:
-                typedef typename std::list<Block>::iterator iterator;
-                typedef typename std::list<Block>::const_iterator const_iterator;
-                
-                // Creates a new block with the given begin and end.
-                Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::size_t id);
-                
-                // Prints the block.
-                void print(Partition const& partition) const;
-
-                // Sets the beginning index of the block.
-                void setBegin(storm::storage::sparse::state_type begin);
-
-                // Moves the beginning index of the block one step further.
-                void incrementBegin();
-                
-                // Sets the end index of the block.
-                void setEnd(storm::storage::sparse::state_type end);
-
-                // Returns the beginning index of the block.
-                storm::storage::sparse::state_type getBegin() const;
-                
-                // Returns the beginning index of the block.
-                storm::storage::sparse::state_type getEnd() const;
-                
-                // Retrieves the original beginning index of the block in case the begin index has been moved.
-                storm::storage::sparse::state_type getOriginalBegin() const;
-                
-                // Returns the iterator the block in the list of overall blocks.
-                iterator getIterator() const;
-
-                // Returns the iterator the block in the list of overall blocks.
-                void setIterator(iterator it);
-                
-                // Returns the iterator the next block in the list of overall blocks if it exists.
-                iterator getNextIterator() const;
-
-                // Returns the iterator the next block in the list of overall blocks if it exists.
-                iterator getPreviousIterator() const;
-
-                // Gets the next block (if there is one).
-                Block& getNextBlock();
-
-                // Gets the next block (if there is one).
-                Block const& getNextBlock() const;
-                
-                // Gets a pointer to the next block (if there is one).
-                Block* getNextBlockPointer();
-
-                // Retrieves whether the block as a successor block.
-                bool hasNextBlock() const;
-
-                // Gets the previous block (if there is one).
-                Block& getPreviousBlock();
-
-                // Gets a pointer to the previous block (if there is one).
-                Block* getPreviousBlockPointer();
-                
-                // Gets the next block (if there is one).
-                Block const& getPreviousBlock() const;
-
-                // Retrieves whether the block as a successor block.
-                bool hasPreviousBlock() const;
-                
-                // Checks consistency of the information in the block.
-                bool check() const;
-                
-                // Retrieves the number of states in this block.
-                std::size_t getNumberOfStates() const;
-                
-                // Checks whether the block is marked as a splitter.
-                bool isMarkedAsSplitter() const;
-                
-                // Marks the block as being a splitter.
-                void markAsSplitter();
-                
-                // Removes the mark.
-                void unmarkAsSplitter();
-                
-                // Retrieves the ID of the block.
-                std::size_t getId() const;
-                
-                // Retrieves the marked position in the block.
-                storm::storage::sparse::state_type getMarkedPosition() const;
-
-                // Sets the marked position to the given value..
-                void setMarkedPosition(storm::storage::sparse::state_type position);
-
-                // Increases the marked position by one.
-                void incrementMarkedPosition();
-                
-                // Resets the marked position to the begin of the block.
-                void resetMarkedPosition();
-                
-                // Retrieves whether the block is marked as a predecessor.
-                bool isMarkedAsPredecessor() const;
-                
-                // Marks the block as being a predecessor block.
-                void markAsPredecessorBlock();
-                
-                // Removes the marking.
-                void unmarkAsPredecessorBlock();
-                
-                // Sets whether or not the block is to be interpreted as absorbing.
-                void setAbsorbing(bool absorbing);
-                
-                // Retrieves whether the block is to be interpreted as absorbing.
-                bool isAbsorbing() const;
-                
-                // Sets the representative state of this block
-                void setRepresentativeState(storm::storage::sparse::state_type representativeState);
-                
-                // Retrieves the representative state for this block.
-                bool hasRepresentativeState() const;
-
-                // Retrieves the representative state for this block.
-                storm::storage::sparse::state_type getRepresentativeState() const;
-                
-            private:
-                // An iterator to itself. This is needed to conveniently insert elements in the overall list of blocks
-                // kept by the partition.
-                iterator selfIt;
-                
-                // Pointers to the next and previous block.
-                Block* next;
-                Block* prev;
-                
-                // The begin and end indices of the block in terms of the state vector of the partition.
-                storm::storage::sparse::state_type begin;
-                storm::storage::sparse::state_type end;
-                
-                // A field that can be used for marking the block.
-                bool markedAsSplitter;
-                
-                // A field that can be used for marking the block as a predecessor block.
-                bool markedAsPredecessorBlock;
-                
-                // A position that can be used to store a certain position within the block.
-                storm::storage::sparse::state_type markedPosition;
-                
-                // A flag indicating whether the block is to be interpreted as absorbing or not.
-                bool absorbing;
-                
-                // The ID of the block. This is only used for debugging purposes.
-                std::size_t id;
-                
-                // An optional representative state for the block. If this is set, this state is used to derive the
-                // atomic propositions of the meta state in the quotient model.
-                boost::optional<storm::storage::sparse::state_type> representativeState;
-            };
-            
-            class Partition {
-            public:
-                friend class Block;
-                
-                /*!
-                 * Creates a partition with one block consisting of all the states.
-                 *
-                 * @param numberOfStates The number of states the partition holds.
-                 * @param keepSilentProbabilities A flag indicating whether or not silent probabilities are to be tracked.
-                 */
-                Partition(std::size_t numberOfStates, bool keepSilentProbabilities = false);
-
-                /*!
-                 * Creates a partition with three blocks: one with all phi states, one with all psi states and one with
-                 * all other states. The former two blocks are marked as being absorbing, because their outgoing
-                 * transitions shall not be taken into account for future refinement.
-                 *
-                 * @param numberOfStates The number of states the partition holds.
-                 * @param prob0States The states which have probability 0 of satisfying phi until psi.
-                 * @param prob1States The states which have probability 1 of satisfying phi until psi.
-                 * @param representativeProb1State If the set of prob1States is non-empty, this needs to be a state
-                 * that is representative for this block in the sense that the state representing this block in the quotient
-                 * model will receive exactly the atomic propositions of the representative state.
-                 * @param keepSilentProbabilities A flag indicating whether or not silent probabilities are to be tracked.
-                 */
-                Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, boost::optional<storm::storage::sparse::state_type> representativeProb1State, bool keepSilentProbabilities = false);
-                
-                Partition() = default;
-                Partition(Partition const& other) = default;
-                Partition& operator=(Partition const& other) = default;
-
-                // Define move-construct and move-assignment explicitly to make sure they exist (as they are vital to
-                // keep validity of pointers.
-				Partition(Partition&& other) : blocks(std::move(other.blocks)), numberOfBlocks(other.numberOfBlocks), stateToBlockMapping(std::move(other.stateToBlockMapping)), statesAndValues(std::move(other.statesAndValues)), positions(std::move(other.positions)), keepSilentProbabilities(other.keepSilentProbabilities), silentProbabilities(std::move(other.silentProbabilities)) {
-                    // Intentionally left empty.
-				}
-                
-				Partition& operator=(Partition&& other) {
-                    if (this != &other) {
-                        blocks = std::move(other.blocks);
-                        numberOfBlocks = other.numberOfBlocks;
-                        stateToBlockMapping = std::move(other.stateToBlockMapping);
-                        statesAndValues = std::move(other.statesAndValues);
-                        positions = std::move(other.positions);
-                        keepSilentProbabilities = other.keepSilentProbabilities;
-                        silentProbabilities = std::move(other.silentProbabilities);
-                    }
-
-					return *this;
-				}
-                
-                /*!
-                 * Splits all blocks of the partition such that afterwards all blocks contain only states with the label
-                 * or no labeled state at all.
-                 */
-                void splitLabel(storm::storage::BitVector const& statesWithLabel);
-                
-                // Retrieves the size of the partition, i.e. the number of blocks.
-                std::size_t size() const;
-                
-                // Prints the partition to the standard output.
-                void print() const;
-
-                // Splits the block at the given position and inserts a new block after the current one holding the rest
-                // of the states.
-                Block& splitBlock(Block& block, storm::storage::sparse::state_type position);
-                
-                // Inserts a block before the given block. The new block will cover all states between the beginning
-                // of the given block and the end of the previous block.
-                Block& insertBlock(Block& block);
-                
-                // Retrieves the blocks of the partition.
-                std::list<Block> const& getBlocks() const;
-
-                // Retrieves the blocks of the partition.
-                std::list<Block>& getBlocks();
-
-                // Checks the partition for internal consistency.
-                bool check() const;
-                
-                // Returns an iterator to the beginning of the states of the given block.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator getBegin(Block const& block);
-                
-                // Returns an iterator to the beginning of the states of the given block.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator getEnd(Block const& block);
-
-                // Returns an iterator to the beginning of the states of the given block.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator getBegin(Block const& block) const;
-                
-                // Returns an iterator to the beginning of the states of the given block.
-                typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::const_iterator getEnd(Block const& block) const;
-
-                // Swaps the positions of the two given states.
-                void swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2);
-
-                // Swaps the positions of the two states given by their positions.
-                void swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2);
-
-                // Retrieves the block of the given state.
-                Block& getBlock(storm::storage::sparse::state_type state);
-
-                // Retrieves the block of the given state.
-                Block const& getBlock(storm::storage::sparse::state_type state) const;
-
-                // Retrieves the position of the given state.
-                storm::storage::sparse::state_type const& getPosition(storm::storage::sparse::state_type state) const;
-
-                // Retrieves the position of the given state.
-                void setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position);
-                
-                // Sets the position of the state to the given position.
-                storm::storage::sparse::state_type const& getState(storm::storage::sparse::state_type position) const;
-
-                // Retrieves the value for the given state.
-                ValueType const& getValue(storm::storage::sparse::state_type state) const;
-                
-                // Retrieves the value at the given position.
-                ValueType const& getValueAtPosition(storm::storage::sparse::state_type position) const;
-                
-                // Sets the given value for the given state.
-                void setValue(storm::storage::sparse::state_type state, ValueType value);
-                
-                // Retrieves the vector with the probabilities going into the current splitter.
-                std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>& getStatesAndValues();
-
-                // Increases the value for the given state by the specified amount.
-                void increaseValue(storm::storage::sparse::state_type state, ValueType value);
-                
-                // Updates the block mapping for the given range of states to the specified block.
-                void updateBlockMapping(Block& block, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last);
-                
-                // Retrieves the first block of the partition.
-                Block& getFirstBlock();
-                
-                // Retrieves whether the given state is fully silent (only in case the silent probabilities are tracked).
-                bool isSilent(storm::storage::sparse::state_type state, storm::utility::ConstantsComparator<ValueType> const& comparator) const;
-                
-                // Retrieves whether the given state has a non-zero silent probability.
-                bool hasSilentProbability(storm::storage::sparse::state_type state, storm::utility::ConstantsComparator<ValueType> const& comparator) const;
-                
-                // Retrieves the silent probability (i.e. the probability to stay within the own equivalence class).
-                ValueType const& getSilentProbability(storm::storage::sparse::state_type state) const;
-                
-                // Sets the silent probabilities for all the states in the range to their values in the range.
-                void setSilentProbabilities(typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last);
-                
-                // Sets the silent probabilities for all states in the range to zero.
-                void setSilentProbabilitiesToZero(typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator first, typename std::vector<std::pair<storm::storage::sparse::state_type, ValueType>>::iterator last);
-
-                // Sets the silent probability for the given state to the given value.
-                void setSilentProbability(storm::storage::sparse::state_type state, ValueType const& value);
-
-            private:
-                // The list of blocks in the partition.
-                std::list<Block> blocks;
-                
-                // The number of blocks of the partition. Although this could be retrieved via block.size(), we want to
-                // guarantee constant access time (which is currently not guaranteed by gcc'c standard libary).
-                uint_fast64_t numberOfBlocks;
-                
-                // A mapping of states to their blocks.
-                std::vector<Block*> stateToBlockMapping;
-                
-                // A vector containing all the states and their values. It is ordered in a special way such that the
-                // blocks only need to define their start/end indices.
-                std::vector<std::pair<storm::storage::sparse::state_type, ValueType>> statesAndValues;
-                
-                // This vector keeps track of the position of each state in the state vector.
-                std::vector<storm::storage::sparse::state_type> positions;
-                
-                // A flag that indicates whether or not the vector with silent probabilities exists.
-                bool keepSilentProbabilities;
-                
-                // This vector keeps track of the silent probabilities (i.e. the probabilities going into the very own
-                // equivalence class) of each state. This means that a state is silent iff its entry is non-zero.
-                std::vector<ValueType> silentProbabilities;
-            };
-            
-            /*!
-             * Performs the partition refinement on the model and thereby computes the equivalence classes under strong
-             * bisimulation equivalence. If required, the quotient model is built and may be retrieved using
-             * getQuotient().
-             *
-             * @param model The model on whose state space to compute the coarses strong bisimulation relation.
-             * @param atomicPropositions The set of atomic propositions that the bisimulation considers.
-             * @param backwardTransitions The backward transitions of the model.
-             * @param The initial partition.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param buildQuotient If set, the quotient model is built and may be retrieved using the getQuotient()
-             * method.
-             * @param comparator A comparator used for comparing constants.
-             */
-            template<typename ModelType>
-            void partitionRefinement(ModelType const& model, std::set<std::string> const& atomicPropositions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition, BisimulationType bisimulationType, bool keepRewards, bool buildQuotient, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            /*!
-             * Refines the partition based on the provided splitter. After calling this method all blocks are stable
-             * with respect to the splitter.
-             *
-             * @param forwardTransitions The forward transitions of the model.
-             * @param backwardTransitions A matrix that can be used to retrieve the predecessors (and their
-             * probabilities).
-             * @param splitter The splitter to use.
-             * @param partition The partition to split.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param splitterQueue A queue into which all blocks that were split are inserted so they can be treated
-             * as splitters in the future.
-             * @param comparator A comparator used for comparing constants.
-             */
-            void refinePartition(storm::storage::SparseMatrix<ValueType> const& forwardTransitions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Block& splitter, Partition& partition, BisimulationType bisimulationType, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            /*!
-             * Refines the block based on their probability values (leading into the splitter).
-             *
-             * @param block The block to refine.
-             * @param partition The partition that contains the block.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param splitterQueue A queue into which all blocks that were split are inserted so they can be treated
-             * as splitters in the future.
-             * @param comparator A comparator used for comparing constants.
-             */
-            void refineBlockProbabilities(Block& block, Partition& partition, BisimulationType bisimulationType, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            void refineBlockWeak(Block& block, Partition& partition, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::deque<Block*>& splitterQueue, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            /*!
-             * Determines the split offsets in the given block.
-             *
-             * @param block The block that is to be analyzed for splits.
-             * @param partition The partition that contains the block.
-             * @param comparator A comparator used for comparing constants.
-             */
-            std::vector<uint_fast64_t> getSplitPointsWeak(Block& block, Partition& partition, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            /*!
-             * Builds the quotient model based on the previously computed equivalence classes (stored in the blocks
-             * of the decomposition.
-             *
-             * @param model The model whose state space was used for computing the equivalence classes. This is used for
-             * determining the transitions of each equivalence class.
-             * @param selectedAtomicPropositions The set of atomic propositions that was considered by the bisimulation.
-             * The quotient will only have these atomic propositions.
-             * @param partition The previously computed partition. This is used for quickly retrieving the block of a
-             * state.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param comparator A comparator used for comparing constants.
-             */
-            template<typename ModelType>
-            void buildQuotient(ModelType const& model, std::set<std::string> const& selectedAtomicPropositions, Partition const& partition, BisimulationType bisimulationType, bool keepRewards, storm::utility::ConstantsComparator<ValueType> const& comparator);
-
-            /*!
-             * Creates the measure-driven initial partition for reaching psi states from phi states.
-             *
-             * @param model The model whose state space is partitioned based on reachability of psi states from phi
-             * states.
-             * @param backwardTransitions The backward transitions of the model.
-             * @param phiStates The phi states in the model.
-             * @param psiStates The psi states in the model.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param bounded If set to true, the initial partition will be chosen in such a way that preserves bounded
-             * reachability queries.
-             * @param comparator A comparator used for comparing constants.
-             * @return The resulting partition.
-             */
-            template<typename ModelType>
-            Partition getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, BisimulationType bisimulationType, bool keepRewards = true, bool bounded = false, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
-            
-            /*!
-             * Creates the initial partition based on all the labels in the given model.
-             *
-             * @param model The model whose state space is partitioned based on its labels.
-             * @param backwardTransitions The backward transitions of the model.
-             * @param bisimulationType The kind of bisimulation that is to be computed.
-             * @param atomicPropositions The set of atomic propositions to respect. If not given, then all atomic
-             * propositions of the model are respected.
-             * @param comparator A comparator used for comparing constants.
-             * @return The resulting partition.
-             */
-            template<typename ModelType>
-            Partition getLabelBasedInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, BisimulationType bisimulationType, boost::optional<std::set<std::string>> const& atomicPropositions = boost::optional<std::set<std::string>>(), bool keepRewards = true, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
-            
-            /*!
-             * Splits all blocks of the given partition into a block that contains all divergent states and another block
-             * containing the non-divergent states.
-             *
-             * @param model The model from which to look-up the probabilities.
-             * @param backwardTransitions The backward transitions of the model.
-             * @param partition The partition that holds the silent probabilities.
-             */
-            template<typename ModelType>
-            void splitOffDivergentStates(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition);
-            
-            /*!
-             * Initializes the silent probabilities by traversing all blocks and adding the probability of going to
-             * the very own equivalence class for each state.
-             *
-             * @param model The model from which to look-up the probabilities.
-             * @param partition The partition that holds the silent probabilities.
-             */
-            template<typename ModelType>
-            void initializeSilentProbabilities(ModelType const& model, Partition& partition);
-            
-            /*!
-             * Splits all blocks of the partition in a way such that the states of each block agree on the rewards.
-             *
-             * @param stateRewardVector The state reward vector of the model.
-             * @param partition The partition that is to be split.
-             * @param comparator A comparator used for comparing constants.
-             */
-            void splitRewards(std::vector<ValueType> const& stateRewardVector, Partition& partition, storm::utility::ConstantsComparator<ValueType> const& comparator);
-            
-            // If required, a quotient model is built and stored in this member.
-            std::shared_ptr<storm::models::sparse::DeterministicModel<ValueType>> quotient;
-        };
-    }
-}
-
-#endif /* STORM_STORAGE_DeterministicModelBisimulationDecomposition_H_ */
\ No newline at end of file
diff --git a/src/storage/Distribution.cpp b/src/storage/Distribution.cpp
index ce42b6ab4..50d794d01 100644
--- a/src/storage/Distribution.cpp
+++ b/src/storage/Distribution.cpp
@@ -3,20 +3,26 @@
 #include <algorithm>
 #include <iostream>
 
+#include "src/utility/macros.h"
+#include "src/utility/constants.h"
+#include "src/utility/ConstantsComparator.h"
+
 #include "src/settings/SettingsManager.h"
 
+#include "src/adapters/CarlAdapter.h"
+
 namespace storm {
     namespace storage {
         
         template<typename ValueType>
-        Distribution<ValueType>::Distribution() : hash(0) {
+        Distribution<ValueType>::Distribution() {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        bool Distribution<ValueType>::operator==(Distribution<ValueType> const& other) const {
+        bool Distribution<ValueType>::equals(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
             // We need to check equality by ourselves, because we need to account for epsilon differences.
-            if (this->distribution.size() != other.distribution.size() || this->getHash() != other.getHash()) {
+            if (this->distribution.size() != other.distribution.size()) {
                 return false;
             }
             
@@ -28,7 +34,7 @@ namespace storm {
                 if (first1->first != first2->first) {
                     return false;
                 }
-                if (std::abs(first1->second - first2->second) > 1e-6) {
+                if (!comparator.isEqual(first1->second, first2->second)) {
                     return false;
                 }
             }
@@ -38,10 +44,28 @@ namespace storm {
         
         template<typename ValueType>
         void Distribution<ValueType>::addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability) {
-            if (this->distribution.find(state) == this->distribution.end()) {
-                this->hash += static_cast<std::size_t>(state);
+            auto it = this->distribution.find(state);
+            if (it == this->distribution.end()) {
+                this->distribution.emplace_hint(it, state, probability);
+            } else {
+                it->second += probability;
+            }
+        }
+        
+        template<typename ValueType>
+        void Distribution<ValueType>::removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator) {
+            auto it = this->distribution.find(state);
+            STORM_LOG_ASSERT(it != this->distribution.end(), "Cannot remove probability, because the state is not in the support of the distribution.");
+            it->second -= probability;
+            if (comparator.isZero(it->second)) {
+                this->distribution.erase(it);
             }
-            this->distribution[state] += probability;
+        }
+        
+        template<typename ValueType>
+        void Distribution<ValueType>::shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator) {
+            removeProbability(fromState, probability, comparator);
+            addProbability(toState, probability);
         }
         
         template<typename ValueType>
@@ -68,7 +92,7 @@ namespace storm {
         void Distribution<ValueType>::scale(storm::storage::sparse::state_type const& state) {
             auto probabilityIterator = this->distribution.find(state);
             if (probabilityIterator != this->distribution.end()) {
-                ValueType scaleValue = 1 / probabilityIterator->second;
+                ValueType scaleValue = storm::utility::one<ValueType>() / probabilityIterator->second;
                 this->distribution.erase(probabilityIterator);
                 
                 for (auto& entry : this->distribution) {
@@ -76,10 +100,10 @@ namespace storm {
                 }
             }
         }
-        
+                
         template<typename ValueType>
-        std::size_t Distribution<ValueType>::getHash() const {
-            return this->hash ^ (this->distribution.size() << 8);
+        std::size_t Distribution<ValueType>::size() const {
+            return this->distribution.size();
         }
         
         template<typename ValueType>
@@ -93,7 +117,34 @@ namespace storm {
             return out;
         }
         
+        template<typename ValueType>
+        bool Distribution<ValueType>::less(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
+            if (this->size() != other.size()) {
+                return this->size() < other.size();
+            }
+            auto firstIt = this->begin();
+            auto firstIte = this->end();
+            auto secondIt = other.begin();
+            for (; firstIt != firstIte; ++firstIt, ++secondIt) {
+                // If the two blocks already differ, we can decide which distribution is smaller.
+                if (firstIt->first != secondIt->first) {
+                    return firstIt->first < secondIt->first;
+                }
+                
+                // If the blocks are the same, but the probability differs, we can also decide which distribution is smaller.
+                if (!comparator.isEqual(firstIt->second, secondIt->second)) {
+                    return comparator.isLess(firstIt->second, secondIt->second);
+                }
+            }
+            return false;
+        }
+        
         template class Distribution<double>;
         template std::ostream& operator<<(std::ostream& out, Distribution<double> const& distribution);
+        
+#ifdef STORM_HAVE_CARL
+        template class Distribution<storm::RationalFunction>;
+        template std::ostream& operator<<(std::ostream& out, Distribution<storm::RationalFunction> const& distribution);
+#endif
     }
 }
diff --git a/src/storage/Distribution.h b/src/storage/Distribution.h
index 010cf43bc..9d3276941 100644
--- a/src/storage/Distribution.h
+++ b/src/storage/Distribution.h
@@ -5,9 +5,16 @@
 #include <ostream>
 #include <boost/container/flat_map.hpp>
 
+#include "src/utility/OsDetection.h"
+
 #include "src/storage/sparse/StateType.h"
 
 namespace storm {
+    namespace utility {
+        template <typename ValueType>
+        class ConstantsComparator;
+    }
+    
     namespace storage {
         
         template<typename ValueType>
@@ -21,6 +28,14 @@ namespace storm {
              * Creates an empty distribution.
              */
             Distribution();
+
+            Distribution(Distribution const& other) = default;
+            Distribution& operator=(Distribution const& other) = default;
+
+#ifndef WINDOWS
+            Distribution(Distribution&& other) = default;
+            Distribution& operator=(Distribution&& other) = default;
+#endif
             
             /*!
              * Checks whether the two distributions specify the same probabilities to go to the same states.
@@ -28,7 +43,7 @@ namespace storm {
              * @param other The distribution with which the current distribution is to be compared.
              * @return True iff the two distributions are equal.
              */
-            bool operator==(Distribution const& other) const;
+            bool equals(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>()) const;
             
             /*!
              * Assigns the given state the given probability under this distribution.
@@ -38,13 +53,34 @@ namespace storm {
              */
             void addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability);
             
+            /*!
+             * Removes the given probability mass of going to the given state.
+             *
+             * @param state The state for which to remove the probability.
+             * @param probability The probability to remove.
+             * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the
+             * entry is removed.
+             */
+            void removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
+            
+            /*!
+             * Removes the probability mass from one state and adds it to another.
+             *
+             * @param fromState The state from which to take the probability mass.
+             * @param toState The state from which to which to add the probability mass.
+             * @param probability The probability mass to shift.
+             * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the
+             * entry is removed.
+             */
+            void shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
+            
             /*!
              * Retrieves an iterator to the elements in this distribution.
              *
              * @return The iterator to the elements in this distribution.
              */
             iterator begin();
-
+            
             /*!
              * Retrieves an iterator to the elements in this distribution.
              *
@@ -58,7 +94,7 @@ namespace storm {
              * @return The iterator past the elements in this distribution.
              */
             iterator end();
-
+            
             /*!
              * Retrieves an iterator past the elements in this distribution.
              *
@@ -76,22 +112,19 @@ namespace storm {
             void scale(storm::storage::sparse::state_type const& state);
             
             /*!
-             * Retrieves the hash value of the distribution.
-             *
-             * @return The hash value of the distribution.
+             * Retrieves the size of the distribution, i.e. the size of the support set.
              */
-            std::size_t getHash() const;
+            std::size_t size() const;
+            
+            bool less(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const;
             
         private:
             // A list of states and the probabilities that are assigned to them.
             container_type distribution;
-            
-            // A hash value that is maintained to allow for quicker equality comparison between distribution.s
-            std::size_t hash;
         };
         
         template<typename ValueType>
-            std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution);
+        std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution);
     }
 }
 
diff --git a/src/storage/bisimulation/BisimulationDecomposition.cpp b/src/storage/bisimulation/BisimulationDecomposition.cpp
new file mode 100644
index 000000000..ae9b61e3a
--- /dev/null
+++ b/src/storage/bisimulation/BisimulationDecomposition.cpp
@@ -0,0 +1,316 @@
+#include "src/storage/bisimulation/BisimulationDecomposition.h"
+
+#include <chrono>
+
+#include "src/models/sparse/Dtmc.h"
+#include "src/models/sparse/Ctmc.h"
+#include "src/models/sparse/Mdp.h"
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/GeneralSettings.h"
+
+#include "src/utility/macros.h"
+#include "src/exceptions/IllegalFunctionCallException.h"
+#include "src/exceptions/InvalidOptionException.h"
+
+namespace storm {
+    namespace storage {
+        
+        using namespace bisimulation;
+        
+        template<typename ModelType, typename BlockDataType>
+        BisimulationDecomposition<ModelType, BlockDataType>::Options::Options(ModelType const& model, storm::logic::Formula const& formula) : Options() {
+            this->preserveSingleFormula(model, formula);
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        BisimulationDecomposition<ModelType, BlockDataType>::Options::Options(ModelType const& model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) : Options() {
+            if (formulas.empty()) {
+                this->respectedAtomicPropositions = model.getStateLabeling().getLabels();
+                this->keepRewards = true;
+            } if (formulas.size() == 1) {
+                this->preserveSingleFormula(model, *formulas.front());
+            } else {
+                for (auto const& formula : formulas) {
+                    preserveFormula(model, *formula);
+                }
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        BisimulationDecomposition<ModelType, BlockDataType>::Options::Options() : measureDrivenInitialPartition(false), phiStates(), psiStates(), respectedAtomicPropositions(), keepRewards(false), type(BisimulationType::Strong), bounded(false), buildQuotient(true) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::Options::preserveFormula(ModelType const& model, storm::logic::Formula const& formula) {
+            // Disable the measure driven initial partition.
+            measureDrivenInitialPartition = false;
+            phiStates = boost::none;
+            psiStates = boost::none;
+            
+            // Preserve rewards if necessary.
+            keepRewards = keepRewards || formula.containsRewardOperator();
+            
+            // Preserve bounded properties if necessary.
+            bounded = bounded || (formula.containsBoundedUntilFormula() || formula.containsNextFormula());
+            
+            // Compute the relevant labels and expressions.
+            this->addToRespectedAtomicPropositions(formula.getAtomicExpressionFormulas(), formula.getAtomicLabelFormulas());
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::Options::preserveSingleFormula(ModelType const& model, storm::logic::Formula const& formula) {
+            keepRewards = formula.containsRewardOperator();
+            
+            // We need to preserve bounded properties iff the formula contains a bounded until or a next subformula.
+            bounded = formula.containsBoundedUntilFormula() || formula.containsNextFormula();
+            
+            // Compute the relevant labels and expressions.
+            this->addToRespectedAtomicPropositions(formula.getAtomicExpressionFormulas(), formula.getAtomicLabelFormulas());
+
+            // Check whether measure driven initial partition is possible and, if so, set it.
+            this->checkAndSetMeasureDrivenInitialPartition(model, formula);
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::Options::checkAndSetMeasureDrivenInitialPartition(ModelType const& model, storm::logic::Formula const& formula) {
+            std::shared_ptr<storm::logic::Formula const> newFormula = formula.asSharedPointer();
+            
+            if (formula.isProbabilityOperatorFormula()) {
+                if (formula.asProbabilityOperatorFormula().hasOptimalityType()) {
+                    optimalityType = formula.asProbabilityOperatorFormula().getOptimalityType();
+                } else if (formula.asProbabilityOperatorFormula().hasBound()) {
+                    storm::logic::ComparisonType comparisonType = formula.asProbabilityOperatorFormula().getComparisonType();
+                    if (comparisonType == storm::logic::ComparisonType::Less || comparisonType == storm::logic::ComparisonType::LessEqual) {
+                        optimalityType = OptimizationDirection::Maximize;
+                    } else {
+                        optimalityType = OptimizationDirection::Minimize;
+                    }
+                }
+                newFormula = formula.asProbabilityOperatorFormula().getSubformula().asSharedPointer();
+            } else if (formula.isRewardOperatorFormula()) {
+                optimalityType = formula.asRewardOperatorFormula().getOptimalityType();
+                newFormula = formula.asRewardOperatorFormula().getSubformula().asSharedPointer();
+            }
+            
+            std::shared_ptr<storm::logic::Formula const> leftSubformula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+            std::shared_ptr<storm::logic::Formula const> rightSubformula;
+            if (newFormula->isUntilFormula()) {
+                leftSubformula = newFormula->asUntilFormula().getLeftSubformula().asSharedPointer();
+                rightSubformula = newFormula->asUntilFormula().getRightSubformula().asSharedPointer();
+                if (leftSubformula->isPropositionalFormula() && rightSubformula->isPropositionalFormula()) {
+                    measureDrivenInitialPartition = true;
+                }
+            } else if (newFormula->isEventuallyFormula()) {
+                rightSubformula = newFormula->asEventuallyFormula().getSubformula().asSharedPointer();
+                if (rightSubformula->isPropositionalFormula()) {
+                    measureDrivenInitialPartition = true;
+                }
+            } else if (newFormula->isReachabilityRewardFormula()) {
+                rightSubformula = newFormula->asReachabilityRewardFormula().getSubformula().asSharedPointer();
+                if (rightSubformula->isPropositionalFormula()) {
+                    measureDrivenInitialPartition = true;
+                }
+            }
+            
+            if (measureDrivenInitialPartition) {
+                storm::modelchecker::SparsePropositionalModelChecker<ModelType> checker(model);
+                std::unique_ptr<storm::modelchecker::CheckResult> phiStatesCheckResult = checker.check(*leftSubformula);
+                std::unique_ptr<storm::modelchecker::CheckResult> psiStatesCheckResult = checker.check(*rightSubformula);
+                phiStates = phiStatesCheckResult->asExplicitQualitativeCheckResult().getTruthValuesVector();
+                psiStates = psiStatesCheckResult->asExplicitQualitativeCheckResult().getTruthValuesVector();
+            } else {
+                optimalityType = boost::none;
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::Options::addToRespectedAtomicPropositions(std::vector<std::shared_ptr<storm::logic::AtomicExpressionFormula const>> const& expressions, std::vector<std::shared_ptr<storm::logic::AtomicLabelFormula const>> const& labels) {
+            std::set<std::string> labelsToRespect;
+            for (auto const& labelFormula : labels) {
+                labelsToRespect.insert(labelFormula->getLabel());
+            }
+            for (auto const& expressionFormula : expressions) {
+                labelsToRespect.insert(expressionFormula->toString());
+            }
+            if (!respectedAtomicPropositions) {
+                respectedAtomicPropositions = labelsToRespect;
+            } else {
+                respectedAtomicPropositions.get().insert(labelsToRespect.begin(), labelsToRespect.end());
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        BisimulationDecomposition<ModelType, BlockDataType>::BisimulationDecomposition(ModelType const& model, Options const& options) : BisimulationDecomposition(model, model.getBackwardTransitions(), options) {
+            // Intentionally left empty.
+        }
+
+        template<typename ModelType, typename BlockDataType>
+        BisimulationDecomposition<ModelType, BlockDataType>::BisimulationDecomposition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Options const& options) : model(model), backwardTransitions(backwardTransitions), options(options), partition(), comparator(), quotient(nullptr) {
+            STORM_LOG_THROW(!options.keepRewards || !model.hasRewardModel() || model.hasUniqueRewardModel(), storm::exceptions::IllegalFunctionCallException, "Bisimulation currently only supports models with at most one reward model.");
+            STORM_LOG_THROW(!options.keepRewards || !model.hasRewardModel() || model.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently supported for models with state rewards only. Consider converting the transition rewards to state rewards (via suitable function calls).");
+            STORM_LOG_THROW(options.type != BisimulationType::Weak || !options.bounded, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation cannot preserve bounded properties.");
+            
+            // Fix the respected atomic propositions if they were not explicitly given.
+            if (!this->options.respectedAtomicPropositions) {
+                this->options.respectedAtomicPropositions = model.getStateLabeling().getLabels();
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::computeBisimulationDecomposition() {
+            std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
+
+            std::chrono::high_resolution_clock::time_point initialPartitionStart = std::chrono::high_resolution_clock::now();
+            // initialize the initial partition.
+            if (options.measureDrivenInitialPartition) {
+                STORM_LOG_THROW(options.phiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without phi states.");
+                STORM_LOG_THROW(options.psiStates, storm::exceptions::InvalidOptionException, "Unable to compute measure-driven initial partition without psi states.");
+                this->initializeMeasureDrivenPartition();
+            } else {
+                this->initializeLabelBasedPartition();
+            }
+            std::chrono::high_resolution_clock::duration initialPartitionTime = std::chrono::high_resolution_clock::now() - initialPartitionStart;
+            
+            this->initialize();
+            
+            std::chrono::high_resolution_clock::time_point refinementStart = std::chrono::high_resolution_clock::now();
+            this->performPartitionRefinement();
+            std::chrono::high_resolution_clock::duration refinementTime = std::chrono::high_resolution_clock::now() - refinementStart;
+            
+            std::chrono::high_resolution_clock::time_point extractionStart = std::chrono::high_resolution_clock::now();
+            this->extractDecompositionBlocks();
+            std::chrono::high_resolution_clock::duration extractionTime = std::chrono::high_resolution_clock::now() - extractionStart;
+            
+            std::chrono::high_resolution_clock::time_point quotientBuildStart = std::chrono::high_resolution_clock::now();
+            if (options.buildQuotient) {
+                this->buildQuotient();
+            }
+            std::chrono::high_resolution_clock::duration quotientBuildTime = std::chrono::high_resolution_clock::now() - quotientBuildStart;
+            
+            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+            
+            if (storm::settings::generalSettings().isShowStatisticsSet()) {
+                std::chrono::milliseconds initialPartitionTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(initialPartitionTime);
+                std::chrono::milliseconds refinementTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(refinementTime);
+                std::chrono::milliseconds extractionTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(extractionTime);
+                std::chrono::milliseconds quotientBuildTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(quotientBuildTime);
+                std::chrono::milliseconds totalTimeInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(totalTime);
+                std::cout << std::endl;
+                std::cout << "Time breakdown:" << std::endl;
+                std::cout << "    * time for initial partition: " << initialPartitionTimeInMilliseconds.count() << "ms" << std::endl;
+                std::cout << "    * time for partitioning: " << refinementTimeInMilliseconds.count() << "ms" << std::endl;
+                std::cout << "    * time for extraction: " << extractionTimeInMilliseconds.count() << "ms" << std::endl;
+                std::cout << "    * time for building quotient: " << quotientBuildTimeInMilliseconds.count() << "ms" << std::endl;
+                std::cout << "------------------------------------------" << std::endl;
+                std::cout << "    * total time: " << totalTimeInMilliseconds.count() << "ms" << std::endl;
+                std::cout << std::endl;
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::performPartitionRefinement() {
+            // Insert all blocks into the splitter queue as a (potential) splitter.
+            std::deque<Block<BlockDataType>*> splitterQueue;
+            std::for_each(partition.getBlocks().begin(), partition.getBlocks().end(), [&] (std::unique_ptr<Block<BlockDataType>> const& block) { block->data().setSplitter(); splitterQueue.push_back(block.get()); } );
+            
+            // Then perform the actual splitting until there are no more splitters.
+            uint_fast64_t iterations = 0;
+            while (!splitterQueue.empty()) {
+                ++iterations;
+
+                // Get and prepare the next splitter.
+                Block<BlockDataType>* splitter = splitterQueue.front();
+                splitterQueue.pop_front();
+                splitter->data().setSplitter(false);
+                
+                // Now refine the partition using the current splitter.
+                refinePartitionBasedOnSplitter(*splitter, splitterQueue);
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        std::shared_ptr<ModelType> BisimulationDecomposition<ModelType, BlockDataType>::getQuotient() const {
+            STORM_LOG_THROW(this->quotient != nullptr, storm::exceptions::IllegalFunctionCallException, "Unable to retrieve quotient model from bisimulation decomposition, because it was not built.");
+            return this->quotient;
+        }
+
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::splitInitialPartitionBasedOnStateRewards() {
+            std::vector<ValueType> const& stateRewardVector = model.getUniqueRewardModel()->second.getStateRewardVector();
+            partition.split([&stateRewardVector] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return stateRewardVector[a] < stateRewardVector[b]; });
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::initializeLabelBasedPartition() {
+            partition = storm::storage::bisimulation::Partition<BlockDataType>(model.getNumberOfStates());
+
+            for (auto const& label : options.respectedAtomicPropositions.get()) {
+                if (label == "init") {
+                    continue;
+                }
+                partition.splitStates(model.getStates(label));
+            }
+            
+            // If the model has state rewards, we need to consider them, because otherwise reward properties are not
+            // preserved.
+            if (options.keepRewards && model.hasRewardModel()) {
+                this->splitInitialPartitionBasedOnStateRewards();
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::initializeMeasureDrivenPartition() {
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = this->getStatesWithProbability01();
+            
+            boost::optional<storm::storage::sparse::state_type> representativePsiState;
+            if (!options.psiStates.get().empty()) {
+                representativePsiState = *options.psiStates.get().begin();
+            }
+            
+            partition = storm::storage::bisimulation::Partition<BlockDataType>(model.getNumberOfStates(), statesWithProbability01.first, options.bounded || options.keepRewards ? options.psiStates.get() : statesWithProbability01.second, representativePsiState);
+            
+            // If the model has state rewards, we need to consider them, because otherwise reward properties are not
+            // preserved.
+            if (options.keepRewards && model.hasRewardModel()) {
+                this->splitInitialPartitionBasedOnStateRewards();
+            }
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::initialize() {
+            // Intentionally left empty.
+        }
+        
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::extractDecompositionBlocks() {
+            // Now move the states from the internal partition into their final place in the decomposition. We do so in
+            // a way that maintains the block IDs as indices.
+            this->blocks.resize(partition.size());
+            for (auto const& blockPtr : partition.getBlocks()) {
+                // We need to sort the states to allow for rapid construction of the blocks.
+                partition.sortBlock(*blockPtr);
+                
+                // Convert the state-value-pairs to states only.
+                this->blocks[blockPtr->getId()] = block_type(partition.begin(*blockPtr), partition.end(*blockPtr), true);
+            }
+        }
+        
+        template class BisimulationDecomposition<storm::models::sparse::Dtmc<double>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Ctmc<double>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Mdp<double>, bisimulation::DeterministicBlockData>;
+        
+#ifdef STORM_HAVE_CARL
+        template class BisimulationDecomposition<storm::models::sparse::Dtmc<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Ctmc<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Mdp<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
+#endif
+    }
+}
diff --git a/src/storage/bisimulation/BisimulationDecomposition.h b/src/storage/bisimulation/BisimulationDecomposition.h
new file mode 100644
index 000000000..a7bd575b2
--- /dev/null
+++ b/src/storage/bisimulation/BisimulationDecomposition.h
@@ -0,0 +1,233 @@
+#ifndef STORM_STORAGE_BISIMULATIONDECOMPOSITION_H_
+#define STORM_STORAGE_BISIMULATIONDECOMPOSITION_H_
+
+#include <deque>
+
+#include "src/storage/sparse/StateType.h"
+#include "src/storage/Decomposition.h"
+#include "src/storage/StateBlock.h"
+#include "src/storage/bisimulation/Partition.h"
+
+#include "src/logic/Formulas.h"
+
+#include "src/utility/constants.h"
+#include "src/utility/ConstantsComparator.h"
+
+namespace storm {
+    namespace utility {
+        template <typename ValueType> class ConstantsComparator;
+    }
+
+    namespace logic {
+        class Formula;
+    }
+    
+    namespace storage {
+
+        enum class BisimulationType { Strong, Weak };
+        
+        /*!
+         * This class is the superclass of all decompositions of a sparse model into its bisimulation quotient.
+         */
+        template <typename ModelType, typename BlockDataType>
+        class BisimulationDecomposition : public Decomposition<StateBlock> {
+        public:
+            typedef typename ModelType::ValueType ValueType;
+            typedef typename ModelType::RewardModelType RewardModelType;
+            
+            // A class that offers the possibility to customize the bisimulation.
+            struct Options {
+                // Creates an object representing the default values for all options.
+                Options();
+                
+                /*!
+                 * Creates an object representing the options necessary to obtain the quotient while still preserving
+                 * the given formula.
+                 *
+                 * @param The model for which the quotient model shall be computed. This needs to be given in order to
+                 * derive a suitable initial partition.
+                 * @param formula The formula that is to be preserved.
+                 */
+                Options(ModelType const& model, storm::logic::Formula const& formula);
+                
+                /*!
+                 * Creates an object representing the options necessary to obtain the smallest quotient while still
+                 * preserving the given formulas.
+                 *
+                 * @param The model for which the quotient model shall be computed. This needs to be given in order to
+                 * derive a suitable initial partition.
+                 * @param formulas The formulas that need to be preserved.
+                 */
+                Options(ModelType const& model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas);
+                
+                /*!
+                 * Changes the options in a way that the given formula is preserved.
+                 *
+                 * @param model The model for which to preserve the formula.
+                 * @param formula The only formula to check.
+                 */
+                void preserveFormula(ModelType const& model, storm::logic::Formula const& formula);
+                
+                // A flag that indicates whether a measure driven initial partition is to be used. If this flag is set
+                // to true, the two optional pairs phiStatesAndLabel and psiStatesAndLabel must be set. Then, the
+                // measure driven initial partition wrt. to the states phi and psi is taken.
+                bool measureDrivenInitialPartition;
+                boost::optional<storm::storage::BitVector> phiStates;
+                boost::optional<storm::storage::BitVector> psiStates;
+                boost::optional<OptimizationDirection> optimalityType;
+                
+                // An optional set of strings that indicate which of the atomic propositions of the model are to be
+                // respected and which may be ignored. If not given, all atomic propositions of the model are respected.
+                boost::optional<std::set<std::string>> respectedAtomicPropositions;
+                
+                // A flag that indicates whether or not the state-rewards of the model are to be respected (and should
+                // be kept in the quotient model, if one is built).
+                bool keepRewards;
+                
+                // A flag that indicates whether a strong or a weak bisimulation is to be computed.
+                BisimulationType type;
+                
+                // A flag that indicates whether step-bounded properties are to be preserved. This may only be set to tru
+                // when computing strong bisimulation equivalence.
+                bool bounded;
+                
+                // A flag that governs whether the quotient model is actually built or only the decomposition is computed.
+                bool buildQuotient;
+                
+            private:
+                /*!
+                 * Sets the options under the assumption that the given formula is the only one that is to be checked.
+                 *
+                 * @param model The model for which to preserve the formula.
+                 * @param formula The only formula to check.
+                 */
+                void preserveSingleFormula(ModelType const& model, storm::logic::Formula const& formula);
+
+                /*!
+                 * Adds the given expressions and labels to the set of respected atomic propositions.
+                 *
+                 * @param expressions The expressions to respect.
+                 * @param labels The labels to respect.
+                 */
+                void addToRespectedAtomicPropositions(std::vector<std::shared_ptr<storm::logic::AtomicExpressionFormula const>> const& expressions, std::vector<std::shared_ptr<storm::logic::AtomicLabelFormula const>> const& labels);
+                
+                /*
+                 * Checks whether a measure driven partition is possible with the given formula and sets the necessary
+                 * data if this is the case.
+                 *
+                 * @param model The model for which to derive the data.
+                 * @param formula The formula for which to derive the data for the measure driven initial partition (if
+                 * applicable).
+                 */
+                void checkAndSetMeasureDrivenInitialPartition(ModelType const& model, storm::logic::Formula const& formula);
+            };
+            
+            /*!
+             * Decomposes the given model into equivalance classes of a bisimulation.
+             *
+             * @param model The model to decompose.
+             * @param options The options to use during for the decomposition.
+             */
+            BisimulationDecomposition(ModelType const& model, Options const& options);
+            
+            /*!
+             * Retrieves the quotient of the model under the computed bisimulation.
+             *
+             * @return The quotient model.
+             */
+            std::shared_ptr<ModelType> getQuotient() const;
+
+            /*!
+             * Computes the decomposition of the model into bisimulation equivalence classes. If requested, a quotient
+             * model is built.
+             */
+            void computeBisimulationDecomposition();
+            
+        protected:
+            /*!
+             * Decomposes the given model into equivalance classes of a bisimulation.
+             *
+             * @param model The model to decompose.
+             * @param backwardTransition The backward transitions of the model.
+             * @param options The options to use during for the decomposition.
+             */
+            BisimulationDecomposition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Options const& options);
+            
+            /*!
+             * Performs the partition refinement on the model and thereby computes the equivalence classes under strong
+             * bisimulation equivalence. If required, the quotient model is built and may be retrieved using
+             * getQuotient().
+             */
+            void performPartitionRefinement();
+            
+            /*!
+             * Refines the partition by considering the given splitter. All blocks that become potential splitters
+             * because of this refinement, are marked as splitters and inserted into the splitter queue.
+             *
+             * @param splitter The splitter to use.
+             * @param splitterQueue The queue into which to insert the newly discovered potential splitters.
+             */
+            virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) = 0;
+            
+            /*!
+             * Builds the quotient model based on the previously computed equivalence classes (stored in the blocks
+             * of the decomposition.
+             */
+            virtual void buildQuotient() = 0;
+            
+            /*!
+             * Initializes the initial partition based on all respected labels.
+             */
+            virtual void initializeLabelBasedPartition();
+            
+            /*!
+             * Creates the measure-driven initial partition for reaching psi states from phi states.
+             */
+            virtual void initializeMeasureDrivenPartition();
+            
+            /*!
+             * A function that can initialize auxiliary data structures. It is called after initializing the initial partition.
+             */
+            virtual void initialize();
+            
+            /*!
+             * Computes the set of states with probability 0/1 for satisfying phi until psi. This is used for the measure
+             * driven initial partition.
+             *
+             * @return The states with probability 0 and 1.
+             */
+            virtual std::pair<storm::storage::BitVector, storm::storage::BitVector> getStatesWithProbability01() = 0;
+            
+            /*!
+             * Splits the initial partition based on the (unique) state reward vector of the model.
+             */
+            virtual void splitInitialPartitionBasedOnStateRewards();
+            
+            /*!
+             * Constructs the blocks of the decomposition object based on the current partition.
+             */
+            void extractDecompositionBlocks();
+            
+            // The model to decompose.
+            ModelType const& model;
+
+            // The backward transitions of the model.
+            storm::storage::SparseMatrix<ValueType> backwardTransitions;
+            
+            // The options used during construction.
+            Options options;
+            
+            // The current partition (used by partition refinement).
+            storm::storage::bisimulation::Partition<BlockDataType> partition;
+            
+            // A comparator used for comparing the distances of constants.
+            storm::utility::ConstantsComparator<ValueType> comparator;
+            
+            // The quotient, if it was build. Otherwhise a null pointer.
+            std::shared_ptr<ModelType> quotient;
+        };
+    }
+}
+
+
+#endif /* STORM_STORAGE_BISIMULATIONDECOMPOSITION_H_ */
diff --git a/src/storage/bisimulation/Block.cpp b/src/storage/bisimulation/Block.cpp
new file mode 100644
index 000000000..812e3b89b
--- /dev/null
+++ b/src/storage/bisimulation/Block.cpp
@@ -0,0 +1,148 @@
+#include "src/storage/bisimulation/Block.h"
+
+#include <iostream>
+#include <iomanip>
+
+#include "src/storage/bisimulation/Partition.h"
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+#include "src/utility/macros.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+
+            template<typename DataType>
+            Block<DataType>::Block(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, Block* previousBlock, Block* nextBlock, std::size_t id) : nextBlock(nextBlock), previousBlock(previousBlock), beginIndex(beginIndex), endIndex(endIndex), id(id), mData() {
+                if (nextBlock != nullptr) {
+                    nextBlock->previousBlock = this;
+                }
+                if (previousBlock != nullptr) {
+                    previousBlock->nextBlock = this;
+                }
+                data().resetMarkers(*this);
+                STORM_LOG_ASSERT(this->beginIndex < this->endIndex, "Unable to create block of illegal size.");
+            }
+            
+            template<typename DataType>
+            bool Block<DataType>::operator==(Block const& other) const {
+                return this == &other;
+            }
+            
+            template<typename DataType>
+            bool Block<DataType>::operator!=(Block const& other) const {
+                return this != &other;
+            }
+            
+            template<typename DataType>
+            void Block<DataType>::print(Partition<DataType> const& partition) const {
+                std::cout << "block [" << this << "] " << this->id << " from " << this->beginIndex << " to " << this->endIndex << " with data [" << this->data() << "]" << std::endl;
+                std::cout << "states ";
+                for (auto stateIt = partition.begin(*this), stateIte = partition.end(*this); stateIt != stateIte; ++stateIt) {
+                    std::cout << *stateIt << ", ";
+                }
+                std::cout << std::endl;
+            }
+            
+            template<typename DataType>
+            void Block<DataType>::setBeginIndex(storm::storage::sparse::state_type beginIndex) {
+                this->beginIndex = beginIndex;
+                data().resetMarkers(*this);
+                STORM_LOG_ASSERT(beginIndex < endIndex, "Unable to resize block to illegal size.");
+            }
+            
+            template<typename DataType>
+            void Block<DataType>::setEndIndex(storm::storage::sparse::state_type endIndex) {
+                this->endIndex = endIndex;
+                data().resetMarkers(*this);
+                STORM_LOG_ASSERT(beginIndex < endIndex, "Unable to resize block to illegal size.");
+            }
+            
+            template<typename DataType>
+            std::size_t Block<DataType>::getId() const {
+                return this->id;
+            }
+            
+            template<typename DataType>
+            storm::storage::sparse::state_type Block<DataType>::getBeginIndex() const {
+                return this->beginIndex;
+            }
+            
+            template<typename DataType>
+            storm::storage::sparse::state_type Block<DataType>::getEndIndex() const {
+                return this->endIndex;
+            }
+            
+            template<typename DataType>
+            Block<DataType> const& Block<DataType>::getNextBlock() const {
+                return *this->nextBlock;
+            }
+            
+            template<typename DataType>
+            bool Block<DataType>::hasNextBlock() const {
+                return this->nextBlock != nullptr;
+            }
+
+            template<typename DataType>
+            Block<DataType>* Block<DataType>::getNextBlockPointer() {
+                return this->nextBlock;
+            }
+            
+            template<typename DataType>
+            Block<DataType> const* Block<DataType>::getNextBlockPointer() const {
+                return this->nextBlock;
+            }
+
+            template<typename DataType>
+            Block<DataType> const& Block<DataType>::getPreviousBlock() const {
+                return *this->previousBlock;
+            }
+
+            template<typename DataType>
+            Block<DataType>* Block<DataType>::getPreviousBlockPointer() {
+                return this->previousBlock;
+            }
+
+            template<typename DataType>
+            Block<DataType> const* Block<DataType>::getPreviousBlockPointer() const {
+                return this->previousBlock;
+            }
+            
+            template<typename DataType>
+            bool Block<DataType>::hasPreviousBlock() const {
+                return this->previousBlock != nullptr;
+            }
+            
+            template<typename DataType>
+            bool Block<DataType>::check() const {
+                STORM_LOG_ASSERT(this->beginIndex < this->endIndex, "Block has negative size.");
+                STORM_LOG_ASSERT(!this->hasPreviousBlock() || this->getPreviousBlock().getNextBlockPointer() == this, "Illegal previous block.");
+                STORM_LOG_ASSERT(!this->hasNextBlock() || this->getNextBlock().getPreviousBlockPointer() == this, "Illegal next block.");
+                return true;
+            }
+            
+            template<typename DataType>
+            std::size_t Block<DataType>::getNumberOfStates() const {
+                return (this->endIndex - this->beginIndex);
+            }
+                        
+            template<typename DataType>
+            DataType& Block<DataType>::data() {
+                return mData;
+            }
+            
+            template<typename DataType>
+            DataType const& Block<DataType>::data() const {
+                return mData;
+            }
+            
+            template<typename DataType>
+            void Block<DataType>::resetMarkers() {
+                mData.resetMarkers(*this);
+            }
+
+            template class Block<DeterministicBlockData>;
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/storage/bisimulation/Block.h b/src/storage/bisimulation/Block.h
new file mode 100644
index 000000000..f9f46051f
--- /dev/null
+++ b/src/storage/bisimulation/Block.h
@@ -0,0 +1,109 @@
+#ifndef STORM_STORAGE_BISIMULATION_BLOCK_H_
+#define STORM_STORAGE_BISIMULATION_BLOCK_H_
+
+#include <list>
+#include <boost/optional.hpp>
+
+#include "src/storage/sparse/StateType.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+            // Forward-declare partition class.
+            template<typename DataType>
+            class Partition;
+            
+            template<typename DataType>
+            class Block {
+            public:
+                friend class Partition<DataType>;
+                
+                // Creates a new block with the given begin and end.
+                Block(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, Block* previous, Block* next, std::size_t id);
+                                
+                Block() = default;
+                Block(Block const& other) = default;
+                Block& operator=(Block const& other) = default;
+                Block(Block&& other) = default;
+                Block& operator=(Block&& other) = default;
+                
+                bool operator==(Block const& other) const;
+                bool operator!=(Block const& other) const;
+                
+                // Prints the block to the standard output.
+                void print(Partition<DataType> const& partition) const;
+                
+                // Returns the beginning index of the block.
+                storm::storage::sparse::state_type getBeginIndex() const;
+                
+                // Returns the beginning index of the block.
+                storm::storage::sparse::state_type getEndIndex() const;
+                
+                // Gets the next block (if there is one).
+                Block const& getNextBlock() const;
+
+                // Gets a pointer to the next block (if there is one).
+                Block* getNextBlockPointer();
+                
+                // Gets a pointer to the next block (if there is one).
+                Block const* getNextBlockPointer() const;
+                
+                // Retrieves whether the block as a successor block.
+                bool hasNextBlock() const;
+                
+                // Gets the next block (if there is one).
+                Block const& getPreviousBlock() const;
+
+                // Gets a pointer to the previous block (if there is one).
+                Block* getPreviousBlockPointer();
+
+                // Gets a pointer to the previous block (if there is one).
+                Block const* getPreviousBlockPointer() const;
+                
+                // Retrieves whether the block as a successor block.
+                bool hasPreviousBlock() const;
+                
+                // Checks consistency of the information in the block.
+                bool check() const;
+                
+                // Retrieves the number of states in this block.
+                std::size_t getNumberOfStates() const;
+                
+                // Retrieves the additional data associated with this block.
+                DataType& data();
+                
+                // Retrieves the additional data associated with this block.
+                DataType const& data() const;
+                
+                // Resets all markers.
+                void resetMarkers();
+                
+                // Retrieves the ID of the block.
+                std::size_t getId() const;
+                
+            private:
+                // Sets the beginning index of the block.
+                void setBeginIndex(storm::storage::sparse::state_type beginIndex);
+                
+                // Sets the end index of the block.
+                void setEndIndex(storm::storage::sparse::state_type endIndex);
+                
+                // Pointers to the next and previous block.
+                Block* nextBlock;
+                Block* previousBlock;
+                
+                // The begin and end indices of the block in terms of the state vector of the partition.
+                storm::storage::sparse::state_type beginIndex;
+                storm::storage::sparse::state_type endIndex;
+                
+                // The ID of the block. This is only used for debugging purposes.
+                std::size_t id;
+                                
+                // A member that stores additional data that depends on the kind of bisimulation.
+                DataType mData;
+            };
+        }
+    }
+}
+
+#endif /* STORM_STORAGE_BISIMULATION_BLOCK_H_ */
diff --git a/src/storage/bisimulation/DeterministicBlockData.cpp b/src/storage/bisimulation/DeterministicBlockData.cpp
new file mode 100644
index 000000000..3a8edced4
--- /dev/null
+++ b/src/storage/bisimulation/DeterministicBlockData.cpp
@@ -0,0 +1,101 @@
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+#include <iostream>
+
+#include "src/exceptions/InvalidOperationException.h"
+#include "src/utility/macros.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+            
+            DeterministicBlockData::DeterministicBlockData() : DeterministicBlockData(0, 0) {
+                // Intentionally left empty.
+            }
+            
+            DeterministicBlockData::DeterministicBlockData(uint_fast64_t marker1, uint_fast64_t marker2) : valMarker1(marker1), valMarker2(marker2), splitterFlag(false), needsRefinementFlag(false), absorbingFlag(false), valRepresentativeState() {
+                // Intentionally left empty.
+            }
+            
+            uint_fast64_t DeterministicBlockData::marker1() const {
+                return valMarker1;
+            }
+            
+            void DeterministicBlockData::setMarker1(uint_fast64_t newMarker1) {
+                valMarker1 = newMarker1;
+            }
+            
+            void DeterministicBlockData::incrementMarker1() {
+                ++valMarker1;
+            }
+            
+            void DeterministicBlockData::decrementMarker1() {
+                --valMarker1;
+            }
+            
+            uint_fast64_t DeterministicBlockData::marker2() const {
+                return valMarker2;
+            }
+            
+            void DeterministicBlockData::setMarker2(uint_fast64_t newMarker2) {
+                valMarker2 = newMarker2;
+            }
+            
+            void DeterministicBlockData::incrementMarker2() {
+                ++valMarker2;
+            }
+            
+            void DeterministicBlockData::decrementMarker2() {
+                --valMarker2;
+            }
+            
+            bool DeterministicBlockData::resetMarkers(Block<DeterministicBlockData> const& block) {
+                bool result = block.getBeginIndex() != this->valMarker1 || block.getBeginIndex() != this->valMarker2;
+                this->valMarker1 = this->valMarker2 = block.getBeginIndex();
+                return result;
+            }
+            
+            bool DeterministicBlockData::splitter() const {
+                return this->splitterFlag;
+            }
+            
+            void DeterministicBlockData::setSplitter(bool value) {
+                this->splitterFlag = value;
+            }
+                        
+            void DeterministicBlockData::setAbsorbing(bool absorbing) {
+                this->absorbingFlag = absorbing;
+            }
+            
+            bool DeterministicBlockData::absorbing() const {
+                return this->absorbingFlag;
+            }
+            
+            void DeterministicBlockData::setRepresentativeState(storm::storage::sparse::state_type representativeState) {
+                this->valRepresentativeState = representativeState;
+            }
+            
+            bool DeterministicBlockData::hasRepresentativeState() const {
+                return static_cast<bool>(valRepresentativeState);
+            }
+            
+            storm::storage::sparse::state_type DeterministicBlockData::representativeState() const {
+                STORM_LOG_THROW(valRepresentativeState, storm::exceptions::InvalidOperationException, "Unable to retrieve representative state for block.");
+                return valRepresentativeState.get();
+            }
+            
+            bool DeterministicBlockData::needsRefinement() const {
+                return needsRefinementFlag;
+            }
+            
+            void DeterministicBlockData::setNeedsRefinement(bool value) {
+                needsRefinementFlag = value;
+            }
+
+            std::ostream& operator<<(std::ostream& out, DeterministicBlockData const& data) {
+                out << "m1: " << data.marker1() << ", m2: " << data.marker2();
+                return out;
+            }
+        }
+    }
+}
diff --git a/src/storage/bisimulation/DeterministicBlockData.h b/src/storage/bisimulation/DeterministicBlockData.h
new file mode 100644
index 000000000..518004efa
--- /dev/null
+++ b/src/storage/bisimulation/DeterministicBlockData.h
@@ -0,0 +1,88 @@
+#ifndef STORM_STORAGE_BISIMULATION_DETERMINISTICBLOCKDATA_H_
+#define STORM_STORAGE_BISIMULATION_DETERMINISTICBLOCKDATA_H_
+
+#include <cstdint>
+
+#include "src/storage/bisimulation/Block.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+            class DeterministicBlockData {
+            public:
+                DeterministicBlockData();
+                DeterministicBlockData(uint_fast64_t marker1, uint_fast64_t marker2);
+                
+                uint_fast64_t marker1() const;
+                void setMarker1(uint_fast64_t newMarker1);
+                void incrementMarker1();
+                void decrementMarker1();
+
+                uint_fast64_t marker2() const;
+                void setMarker2(uint_fast64_t newMarker2);
+                void incrementMarker2();
+                void decrementMarker2();
+                
+                /*!
+                 * This method needs to be called whenever the block was modified to reset the data of the change.
+                 *
+                 * @param block The block that this data belongs to.
+                 * @return True iff the data changed as a consequence of notifying it.
+                 */
+                bool resetMarkers(Block<DeterministicBlockData> const& block);
+                
+                // Checks whether the block is marked as a splitter.
+                bool splitter() const;
+                
+                // Marks the block as being a splitter.
+                void setSplitter(bool value = true);
+                
+                // Retrieves whether the block is marked as a predecessor.
+                bool needsRefinement() const;
+                
+                // Marks the block as needing refinement (or not).
+                void setNeedsRefinement(bool value = true);
+                
+                // Sets whether or not the block is to be interpreted as absorbing.
+                void setAbsorbing(bool absorbing);
+                
+                // Retrieves whether the block is to be interpreted as absorbing.
+                bool absorbing() const;
+                
+                // Sets the representative state of this block
+                void setRepresentativeState(storm::storage::sparse::state_type representativeState);
+                
+                // Retrieves whether this block has a representative state.
+                bool hasRepresentativeState() const;
+                
+                // Retrieves the representative state for this block.
+                storm::storage::sparse::state_type representativeState() const;
+                
+                friend std::ostream& operator<<(std::ostream& out, DeterministicBlockData const& data);
+                
+            public:
+                // Two markers that can be used for various purposes. Whenever the block is split, both the markers are
+                // set to the beginning index of the block.
+                uint_fast64_t valMarker1;
+                uint_fast64_t valMarker2;
+                
+                // A flag that can be used for marking the block as being a splitter.
+                bool splitterFlag;
+                
+                // A flag that can be used for marking the block as needing refinement.
+                bool needsRefinementFlag;
+                
+                // A flag indicating whether the block is to be interpreted as absorbing or not.
+                bool absorbingFlag;
+                
+                // An optional representative state for the block. If this is set, this state is used to derive the
+                // atomic propositions of the meta state in the quotient model.
+                boost::optional<storm::storage::sparse::state_type> valRepresentativeState;
+            };
+            
+            std::ostream& operator<<(std::ostream& out, DeterministicBlockData const& data);
+        }
+    }
+}
+
+#endif /* STORM_STORAGE_BISIMULATION_DETERMINISTICBLOCKDATA_H_ */
diff --git a/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp
new file mode 100644
index 000000000..f81519084
--- /dev/null
+++ b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp
@@ -0,0 +1,642 @@
+#include "src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h"
+
+#include <algorithm>
+#include <unordered_map>
+#include <chrono>
+#include <iomanip>
+#include <boost/iterator/zip_iterator.hpp>
+
+#include "src/adapters/CarlAdapter.h"
+#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
+
+#include "src/models/sparse/Dtmc.h"
+#include "src/models/sparse/Ctmc.h"
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/utility/graph.h"
+#include "src/utility/constants.h"
+#include "src/utility/ConstantsComparator.h"
+#include "src/exceptions/IllegalFunctionCallException.h"
+#include "src/exceptions/InvalidArgumentException.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/GeneralSettings.h"
+
+namespace storm {
+    namespace storage {
+        
+        using namespace bisimulation;
+        
+        template<typename ModelType>
+        DeterministicModelBisimulationDecomposition<ModelType>::DeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, DeterministicModelBisimulationDecomposition::BlockDataType>::Options const& options) : BisimulationDecomposition<ModelType, DeterministicModelBisimulationDecomposition::BlockDataType>(model, options), probabilitiesToCurrentSplitter(model.getNumberOfStates(), storm::utility::zero<ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ModelType>
+        std::pair<storm::storage::BitVector, storm::storage::BitVector> DeterministicModelBisimulationDecomposition<ModelType>::getStatesWithProbability01() {
+            return storm::utility::graph::performProb01(this->backwardTransitions, this->options.phiStates.get(), this->options.psiStates.get());
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::splitOffDivergentStates() {
+            std::vector<storm::storage::sparse::state_type> stateStack;
+            stateStack.reserve(this->model.getNumberOfStates());
+            storm::storage::BitVector nondivergentStates(this->model.getNumberOfStates());
+            
+            uint_fast64_t currentSize = this->partition.size();
+            for (uint_fast64_t blockIndex = 0; blockIndex < currentSize; ++blockIndex) {
+                auto& block = *this->partition.getBlocks()[blockIndex];
+                nondivergentStates.clear();
+                
+                for (auto stateIt = this->partition.begin(block), stateIte = this->partition.end(block); stateIt != stateIte; ++stateIt) {
+                    if (nondivergentStates.get(*stateIt)) {
+                        continue;
+                    }
+                    
+                    // Now traverse the forward transitions of the current state and check whether there is a
+                    // transition to some other block.
+                    bool isDirectlyNonDivergent = false;
+                    for (auto const& successor : this->model.getRows(*stateIt)) {
+                        // If there is such a transition, then we can mark all states in the current block that can
+                        // reach the state as non-divergent.
+                        if (this->partition.getBlock(successor.getColumn()) != block) {
+                            isDirectlyNonDivergent = true;
+                            break;
+                        }
+                    }
+                    
+                    if (isDirectlyNonDivergent) {
+                        stateStack.push_back(*stateIt);
+                        
+                        while (!stateStack.empty()) {
+                            storm::storage::sparse::state_type currentState = stateStack.back();
+                            stateStack.pop_back();
+                            nondivergentStates.set(currentState);
+                            
+                            for (auto const& predecessor : this->backwardTransitions.getRow(currentState)) {
+                                if (this->partition.getBlock(predecessor.getColumn()) == block && !nondivergentStates.get(predecessor.getColumn())) {
+                                    stateStack.push_back(predecessor.getColumn());
+                                }
+                            }
+                        }
+                    }
+                }
+                
+                if (!nondivergentStates.empty() && nondivergentStates.getNumberOfSetBits() != block.getNumberOfStates()) {
+                    // After performing the split, the current block will contain the divergent states only.
+                    this->partition.splitStates(block, nondivergentStates);
+                    
+                    // Since the remaining states in the block are divergent, we can mark the block as absorbing.
+                    // This also guarantees that the self-loop will be added to the state of the quotient
+                    // representing this block of states.
+                    block.data().setAbsorbing(true);
+                } else if (nondivergentStates.empty()) {
+                    // If there are only diverging states in the block, we need to make it absorbing.
+                    block.data().setAbsorbing(true);
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::initializeSilentProbabilities() {
+            silentProbabilities.resize(this->model.getNumberOfStates(), storm::utility::zero<ValueType>());
+            for (storm::storage::sparse::state_type state = 0; state < this->model.getNumberOfStates(); ++state) {
+                Block<BlockDataType> const* currentBlockPtr = &this->partition.getBlock(state);
+                for (auto const& successorEntry : this->model.getRows(state)) {
+                    if (&this->partition.getBlock(successorEntry.getColumn()) == currentBlockPtr) {
+                        silentProbabilities[state] += successorEntry.getValue();
+                    }
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::initializeWeakDtmcBisimulation() {
+            // If we are creating the initial partition for weak bisimulation on DTMCs, we need to (a) split off all
+            // divergent states of each initial block and (b) initialize the vector of silent probabilities.
+            this->splitOffDivergentStates();
+            this->initializeSilentProbabilities();
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::initializeMeasureDrivenPartition() {
+            BisimulationDecomposition<ModelType, BlockDataType>::initializeMeasureDrivenPartition();
+            
+            if (this->options.type == BisimulationType::Weak && this->model.getType() == storm::models::ModelType::Dtmc) {
+                this->initializeWeakDtmcBisimulation();
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::initializeLabelBasedPartition() {
+            BisimulationDecomposition<ModelType, BlockDataType>::initializeLabelBasedPartition();
+            
+            if (this->options.type == BisimulationType::Weak && this->model.getType() == storm::models::ModelType::Dtmc) {
+                this->initializeWeakDtmcBisimulation();
+            }
+        }
+        
+        template<typename ModelType>
+        typename DeterministicModelBisimulationDecomposition<ModelType>::ValueType const& DeterministicModelBisimulationDecomposition<ModelType>::getProbabilityToSplitter(storm::storage::sparse::state_type const& state) const {
+            return probabilitiesToCurrentSplitter[state];
+        }
+        
+        template<typename ModelType>
+        bool DeterministicModelBisimulationDecomposition<ModelType>::isSilent(storm::storage::sparse::state_type const& state) const {
+            return this->comparator.isOne(silentProbabilities[state]);
+        }
+        
+        template<typename ModelType>
+        bool DeterministicModelBisimulationDecomposition<ModelType>::hasNonZeroSilentProbability(storm::storage::sparse::state_type const& state) const {
+            return !this->comparator.isZero(silentProbabilities[state]);
+        }
+        
+        template<typename ModelType>
+        typename DeterministicModelBisimulationDecomposition<ModelType>::ValueType DeterministicModelBisimulationDecomposition<ModelType>::getSilentProbability(storm::storage::sparse::state_type const& state) const {
+            return silentProbabilities[state];
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitterStrong(std::list<Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            for (auto block : predecessorBlocks) {
+                STORM_LOG_TRACE("Refining predecessor " << block->getId() << " of splitter");
+
+                // Depending on the actions we need to take, the block to refine changes, so we need to keep track of it.
+                Block<BlockDataType>* blockToRefineProbabilistically = block;
+                
+                bool split = false;
+                // If the new begin index has shifted to a non-trivial position, we need to split the block.
+                if (block->getBeginIndex() != block->data().marker1() && block->getEndIndex() != block->data().marker1()) {
+                    split = true;
+                    this->partition.splitBlock(*block, block->data().marker1());
+                    blockToRefineProbabilistically = block->getPreviousBlockPointer();
+                }
+                
+                split |= this->partition.splitBlock(*blockToRefineProbabilistically,
+                                                    [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                                                        return this->comparator.isLess(getProbabilityToSplitter(state1), getProbabilityToSplitter(state2));
+                                                    },
+                                                    [&splitterQueue] (Block<BlockDataType>& block) {
+                                                        splitterQueue.emplace_back(&block); block.data().setSplitter();
+                                                    });
+                
+                
+                // If the predecessor block was split, we need to insert it into the splitter queue if it is not already
+                // marked as a splitter.
+                if (split && !blockToRefineProbabilistically->data().splitter()) {
+                    splitterQueue.emplace_back(blockToRefineProbabilistically);
+                    blockToRefineProbabilistically->data().setSplitter();
+                }
+                
+                // If the block was *not* split, we need to reset the markers by notifying the data.
+                block->resetMarkers();
+                
+                // Remember that we have refined the block.
+                block->data().setNeedsRefinement(false);
+            }
+        }
+        
+        template<typename ModelType>
+        bool DeterministicModelBisimulationDecomposition<ModelType>::possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& predecessorBlock) const {
+            return predecessorBlock.getNumberOfStates() > 1 && !predecessorBlock.data().absorbing();
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::increaseProbabilityToSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType> const& predecessorBlock, ValueType const& value) {
+            storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
+            
+            // If the position of the state is to the right of marker1, we have not seen it before.
+            if (predecessorPosition >= predecessorBlock.data().marker1()) {
+                // Then, we just set the value.
+                probabilitiesToCurrentSplitter[predecessor] = value;
+            } else {
+                // If the state was seen as a predecessor before, we add the value to the existing value.
+                probabilitiesToCurrentSplitter[predecessor] += value;
+            }
+        }
+        
+        template <typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock) {
+            this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1()));
+            predecessorBlock.data().incrementMarker1();
+        }
+        
+        template <typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock) {
+            this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker2()));
+            predecessorBlock.data().incrementMarker2();
+        }
+        
+        template <typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip) {
+            storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
+            
+            // If the predecessors of the given predecessor were already explored, we can move it easily.
+            if (predecessorPosition <= currentPositionInSplitter + elementsToSkip) {
+                this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1()));
+                predecessorBlock.data().incrementMarker1();
+            } else {
+                // Otherwise, we need to move the predecessor, but we need to make sure that we explore its
+                // predecessors later. We do this by moving it to a range at the beginning of the block that will hold
+                // all predecessors in the splitter whose predecessors have yet to be explored.
+                if (predecessorBlock.data().marker2() == predecessorBlock.data().marker1()) {
+                    this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition);
+                    this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1);
+                } else {
+                    this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition);
+                    this->partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.data().marker1());
+                    this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1);
+                }
+                
+                // Since we had to move an already explored state to the right of the current position,
+                ++elementsToSkip;
+                predecessorBlock.data().incrementMarker1();
+                predecessorBlock.data().incrementMarker2();
+            }
+        }
+        
+        template <typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::exploreRemainingStatesOfSplitter(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks) {
+            for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.begin(splitter) + (splitter.data().marker2() - splitter.getBeginIndex()); splitterIt != splitterIte; ++splitterIt) {
+                storm::storage::sparse::state_type currentState = *splitterIt;
+                
+                for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
+                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessor);
+                    
+                    // If the block does not need to be refined, we skip it.
+                    if (!possiblyNeedsRefinement(predecessorBlock)) {
+                        continue;
+                    }
+
+                    // We keep track of the probability of the predecessor moving to the splitter.
+                    increaseProbabilityToSplitter(predecessor, predecessorBlock, predecessorEntry.getValue());
+                    
+                    // Only move the state if it has not been seen as a predecessor before.
+                    storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
+                    if (predecessorPosition >= predecessorBlock.data().marker1()) {
+                        moveStateToMarker1(predecessor, predecessorBlock);
+                    }
+                    
+                    insertIntoPredecessorList(predecessorBlock, predecessorBlocks);
+                }
+            }
+            
+            // Finally, we can reset the second marker.
+            splitter.data().setMarker2(splitter.getBeginIndex());
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(bisimulation::Block<BlockDataType>& block) {
+            // For all predecessors, we can set the probability to the current probability of moving to the splitter.
+            for (auto stateIt = this->partition.begin(block), stateIte = this->partition.begin() + block.data().marker1(); stateIt != stateIte; ++stateIt) {
+                silentProbabilities[*stateIt] = probabilitiesToCurrentSplitter[*stateIt];
+            }
+            // All non-predecessors have a silent probability of zero.
+            for (auto stateIt = this->partition.begin() + block.data().marker1(), stateIte = this->partition.end(block); stateIt != stateIte; ++stateIt) {
+                silentProbabilities[*stateIt] = storm::utility::zero<ValueType>();
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::updateSilentProbabilitiesBasedOnTransitions(bisimulation::Block<BlockDataType>& block) {
+            for (auto stateIt = this->partition.begin(block), stateIte = this->partition.end(block); stateIt != stateIte; ++stateIt) {
+                if (hasNonZeroSilentProbability(*stateIt)) {
+                    ValueType newSilentProbability = storm::utility::zero<ValueType>();
+                    for (auto const& successorEntry : this->model.getTransitionMatrix().getRow(*stateIt)) {
+                        if (this->partition.getBlock(successorEntry.getColumn()) == block) {
+                            newSilentProbability += successorEntry.getValue();
+                        }
+                    }
+                    silentProbabilities[*stateIt] = newSilentProbability;
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::computeConditionalProbabilitiesForNonSilentStates(bisimulation::Block<BlockDataType>& block) {
+            for (auto stateIt = this->partition.begin() + block.getBeginIndex(), stateIte = this->partition.begin() + block.data().marker1(); stateIt != stateIte; ++stateIt) {
+                if (!this->comparator.isOne(getSilentProbability(*stateIt))) {
+                    probabilitiesToCurrentSplitter[*stateIt] /= storm::utility::one<ValueType>() - getSilentProbability(*stateIt);
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        std::vector<uint_fast64_t> DeterministicModelBisimulationDecomposition<ModelType>::computeNonSilentBlocks(bisimulation::Block<BlockDataType>& block) {
+            auto less = [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { return probabilitiesToCurrentSplitter[state1] < probabilitiesToCurrentSplitter[state2]; };
+            this->partition.sortRange(block.getBeginIndex(), block.data().marker1(), less);
+            return this->partition.computeRangesOfEqualValue(block.getBeginIndex(), block.data().marker1(), less);
+        }
+        
+        template<typename ModelType>
+        std::vector<storm::storage::BitVector> DeterministicModelBisimulationDecomposition<ModelType>::computeWeakStateLabelingBasedOnNonSilentBlocks(bisimulation::Block<BlockDataType> const& block, std::vector<uint_fast64_t> const& nonSilentBlockIndices) {
+            // Now that we have the split points of the non-silent states, we perform a backward search from
+            // each non-silent state and label the predecessors with the class of the non-silent state.
+            std::vector<storm::storage::BitVector> stateLabels(block.getNumberOfStates(), storm::storage::BitVector(nonSilentBlockIndices.size() - 1));
+            
+            std::vector<storm::storage::sparse::state_type> stateStack;
+            stateStack.reserve(block.getNumberOfStates());
+            for (uint_fast64_t stateClassIndex = 0; stateClassIndex < nonSilentBlockIndices.size() - 1; ++stateClassIndex) {
+                for (auto stateIt = this->partition.begin() + nonSilentBlockIndices[stateClassIndex], stateIte = this->partition.begin() + nonSilentBlockIndices[stateClassIndex + 1]; stateIt != stateIte; ++stateIt) {
+                    
+                    stateStack.push_back(*stateIt);
+                    stateLabels[this->partition.getPosition(*stateIt) - block.getBeginIndex()].set(stateClassIndex);
+                    while (!stateStack.empty()) {
+                        storm::storage::sparse::state_type currentState = stateStack.back();
+                        stateStack.pop_back();
+                        
+                        for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
+                            storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+
+                            if (this->comparator.isZero(predecessorEntry.getValue())) {
+                                continue;
+                            }
+
+                            // Only if the state is in the same block, is a silent state and it has not yet been
+                            // labeled with the current label.
+                            if (this->partition.getBlock(predecessor) == block && isSilent(predecessor) && !stateLabels[this->partition.getPosition(predecessor) - block.getBeginIndex()].get(stateClassIndex)) {
+                                stateStack.push_back(predecessor);
+                                stateLabels[this->partition.getPosition(predecessor) - block.getBeginIndex()].set(stateClassIndex);
+                            }
+                        }
+                    }
+                }
+            }
+            
+            return stateLabels;
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlockOfSplitterWeak(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            // First, we need to turn the one-step probabilities to go to the splitter to the conditional probabilities
+            // for all non-silent states.
+            computeConditionalProbabilitiesForNonSilentStates(block);
+            
+            // Then, we need to compute a labeling of the states that expresses which of the non-silent blocks they can reach.
+            std::vector<uint_fast64_t> nonSilentBlockIndices = computeNonSilentBlocks(block);
+            std::vector<storm::storage::BitVector> weakStateLabels = computeWeakStateLabelingBasedOnNonSilentBlocks(block, nonSilentBlockIndices);
+            
+            // Then split the block according to this labeling.
+            // CAUTION: that this assumes that the positions of the states in the partition are not update until after
+            // the sorting is over. Otherwise, this interferes with the data used in the sorting process.
+            storm::storage::sparse::state_type originalBlockIndex = block.getBeginIndex();
+            auto split = this->partition.splitBlock(block,
+                                                     [&weakStateLabels,&block,originalBlockIndex,this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                                                         return weakStateLabels[this->partition.getPosition(state1) - originalBlockIndex] < weakStateLabels[this->partition.getPosition(state2) - originalBlockIndex];
+                                                     },
+                                                     [this, &splitterQueue] (bisimulation::Block<BlockDataType>& block) {
+                                                         updateSilentProbabilitiesBasedOnTransitions(block);
+                                                         
+                                                         // Insert the new block as a splitter.
+                                                         block.data().setSplitter();
+                                                         splitterQueue.emplace_back(&block);
+                                                     });
+            
+            // If the block was split, we also update the silent probabilities.
+            if (split) {
+                updateSilentProbabilitiesBasedOnTransitions(block);
+                
+                if (!block.data().splitter()) {
+                    // Insert the new block as a splitter.
+                    block.data().setSplitter();
+                    splitterQueue.emplace_back(&block);
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitterWeak(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            for (auto block : predecessorBlocks) {
+                if (*block != splitter) {
+                    refinePredecessorBlockOfSplitterWeak(*block, splitterQueue);
+                } else {
+                    // If the block to split is the splitter itself, we must not do any splitting here.
+                }
+                
+                block->resetMarkers();
+                block->data().setNeedsRefinement(false);
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::insertIntoPredecessorList(bisimulation::Block<BlockDataType>& predecessorBlock, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks) {
+            // Insert the block into the list of blocks to refine (if that has not already happened).
+            if (!predecessorBlock.data().needsRefinement()) {
+                predecessorBlocks.emplace_back(&predecessorBlock);
+                predecessorBlock.data().setNeedsRefinement();
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            STORM_LOG_TRACE("Refining partition based on splitter " << splitter.getId());
+            
+            // The outline of the refinement is as follows.
+            //
+            // We iterate over all states of the splitter and determine for each predecessor the state the probability
+            // entering the splitter. These probabilities are written to a member vector so that after the iteration
+            // process we have the probabilities of all predecessors of the splitter of entering the splitter in one
+            // step. To directly separate the states having a transition into the splitter from the ones who do not,
+            // we move the states to certain locations. That is, on encountering a predecessor of the splitter, it is
+            // moved to the beginning of its block. If the predecessor is in the splitter itself, we have to be a bit
+            // careful about where to move states.
+            //
+            // After this iteration, there may be states of the splitter whose predecessors have not yet been explored,
+            // so this needs to be done now.
+            //
+            // Finally, we use the information obtained in the first part for the actual splitting process in which all
+            // predecessor blocks of the splitter are split based on the probabilities computed earlier.
+            std::list<Block<BlockDataType>*> predecessorBlocks;
+            storm::storage::sparse::state_type currentPosition = splitter.getBeginIndex();
+            bool splitterIsPredecessorBlock = false;
+            for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.end(splitter); splitterIt != splitterIte; ++splitterIt, ++currentPosition) {
+                storm::storage::sparse::state_type currentState = *splitterIt;
+                
+                uint_fast64_t elementsToSkip = 0;
+                for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
+                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+                    storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
+                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessor);
+                    
+                    // If the block does not need to be refined, we skip it.
+                    if (!possiblyNeedsRefinement(predecessorBlock)) {
+                        continue;
+                    }
+                    
+                    // We keep track of the probability of the predecessor moving to the splitter.
+                    increaseProbabilityToSplitter(predecessor, predecessorBlock, predecessorEntry.getValue());
+                    
+                    // We only need to move the predecessor if its not already known as a predecessor already.
+                    if (predecessorPosition >= predecessorBlock.data().marker1()) {
+                        // If the predecessor block is not the splitter, we can move the state easily.
+                        if (predecessorBlock != splitter) {
+                            moveStateToMarker1(predecessor, predecessorBlock);
+                        } else {
+                            // If the predecessor is in the splitter, we need to be a bit more careful.
+                            splitterIsPredecessorBlock = true;
+                            moveStateInSplitter(predecessor, predecessorBlock, currentPosition, elementsToSkip);
+                        }
+                        
+                        insertIntoPredecessorList(predecessorBlock, predecessorBlocks);
+                    }
+                }
+                
+                // If, as a consequence of shifting states, we need to skip some elements, do so now.
+                splitterIt += elementsToSkip;
+                currentPosition += elementsToSkip;
+            }
+            
+            // If the splitter was a predecessor block of itself, we potentially need to explore some states that have
+            // not been explored yet.
+            if (splitterIsPredecessorBlock) {
+                exploreRemainingStatesOfSplitter(splitter, predecessorBlocks);
+            }
+            
+            // Finally, we split the block based on the precomputed probabilities and the chosen bisimulation type.
+            if (this->options.type == BisimulationType::Strong || this->model.getType() == storm::models::ModelType::Ctmc) {
+                refinePredecessorBlocksOfSplitterStrong(predecessorBlocks, splitterQueue);
+            } else {
+                // If the splitter is a predecessor of we can use the computed probabilities to update the silent
+                // probabilities.
+                if (splitterIsPredecessorBlock) {
+                    updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(splitter);
+                }
+                
+                refinePredecessorBlocksOfSplitterWeak(splitter, predecessorBlocks, splitterQueue);
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
+            // In order to create the quotient model, we need to construct
+            // (a) the new transition matrix,
+            // (b) the new labeling,
+            // (c) the new reward structures.
+            
+            // Prepare a matrix builder for (a).
+            storm::storage::SparseMatrixBuilder<ValueType> builder(this->size(), this->size());
+            
+            // Prepare the new state labeling for (b).
+            storm::models::sparse::StateLabeling newLabeling(this->size());
+            std::set<std::string> atomicPropositionsSet = this->options.respectedAtomicPropositions.get();
+            atomicPropositionsSet.insert("init");
+            std::vector<std::string> atomicPropositions = std::vector<std::string>(atomicPropositionsSet.begin(), atomicPropositionsSet.end());
+            for (auto const& ap : atomicPropositions) {
+                newLabeling.addLabel(ap);
+            }
+            
+            // If the model had state rewards, we need to build the state rewards for the quotient as well.
+            boost::optional<std::vector<ValueType>> stateRewards;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                stateRewards = std::vector<ValueType>(this->blocks.size());
+            }
+            
+            // Now build (a) and (b) by traversing all blocks.
+            for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) {
+                auto const& block = this->blocks[blockIndex];
+                
+                // Pick one representative state. For strong bisimulation it doesn't matter which state it is, because
+                // they all behave equally.
+                storm::storage::sparse::state_type representativeState = *block.begin();
+                
+                // However, for weak bisimulation, we need to make sure the representative state is a non-silent one (if
+                // there is any such state).
+                if (this->options.type == BisimulationType::Weak) {
+                    for (auto const& state : block) {
+                        if (!isSilent(state)) {
+                            representativeState = state;
+                            break;
+                        }
+                    }
+                }
+                
+                Block<BlockDataType> const& oldBlock = this->partition.getBlock(representativeState);
+                
+                // If the block is absorbing, we simply add a self-loop.
+                if (oldBlock.data().absorbing()) {
+                    builder.addNextValue(blockIndex, blockIndex, storm::utility::one<ValueType>());
+                    
+                    // If the block has a special representative state, we retrieve it now.
+                    if (oldBlock.data().hasRepresentativeState()) {
+                        representativeState = oldBlock.data().representativeState();
+                    }
+                    
+                    // Add all of the selected atomic propositions that hold in the representative state to the state
+                    // representing the block.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                } else {
+                    // Compute the outgoing transitions of the block.
+                    std::map<storm::storage::sparse::state_type, ValueType> blockProbability;
+                    for (auto const& entry : this->model.getTransitionMatrix().getRow(representativeState)) {
+                        storm::storage::sparse::state_type targetBlock = this->partition.getBlock(entry.getColumn()).getId();
+                        
+                        // If we are computing a weak bisimulation quotient, there is no need to add self-loops.
+                        if ((this->options.type == BisimulationType::Weak) && targetBlock == blockIndex) {
+                            continue;
+                        }
+                        
+                        auto probIterator = blockProbability.find(targetBlock);
+                        if (probIterator != blockProbability.end()) {
+                            probIterator->second += entry.getValue();
+                        } else {
+                            blockProbability[targetBlock] = entry.getValue();
+                        }
+                    }
+                    
+                    // Now add them to the actual matrix.
+                    for (auto const& probabilityEntry : blockProbability) {
+                        if (this->options.type == BisimulationType::Weak && this->model.getType() == storm::models::ModelType::Dtmc) {
+                            builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second / (storm::utility::one<ValueType>() - getSilentProbability(representativeState)));
+                        } else {
+                            builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second);
+                        }
+                    }
+                    
+                    // Otherwise add all atomic propositions to the equivalence class that the representative state
+                    // satisfies.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                }
+                
+                // If the model has state rewards, we simply copy the state reward of the representative state, because
+                // all states in a block are guaranteed to have the same state reward.
+                if (this->options.keepRewards && this->model.hasRewardModel()) {
+                    typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                    stateRewards.get()[blockIndex] = nameRewardModelPair->second.getStateRewardVector()[representativeState];
+                }
+            }
+            
+            // Now check which of the blocks of the partition contain at least one initial state.
+            for (auto initialState : this->model.getInitialStates()) {
+                Block<BlockDataType> const& initialBlock = this->partition.getBlock(initialState);
+                newLabeling.addLabelToState("init", initialBlock.getId());
+            }
+            
+            // Construct the reward model mapping.
+            std::unordered_map<std::string, typename ModelType::RewardModelType> rewardModels;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                rewardModels.insert(std::make_pair(nameRewardModelPair->first, typename ModelType::RewardModelType(stateRewards)));
+            }
+            
+            // Finally construct the quotient model.
+            this->quotient = std::shared_ptr<ModelType>(new ModelType(builder.build(), std::move(newLabeling), std::move(rewardModels)));
+        }
+        
+        template class DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>;
+        template class DeterministicModelBisimulationDecomposition<storm::models::sparse::Ctmc<double>>;
+        
+#ifdef STORM_HAVE_CARL
+        template class DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<storm::RationalFunction>>;
+        template class DeterministicModelBisimulationDecomposition<storm::models::sparse::Ctmc<storm::RationalFunction>>;
+#endif
+    }
+}
diff --git a/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h
new file mode 100644
index 000000000..b9f2a2a5a
--- /dev/null
+++ b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h
@@ -0,0 +1,130 @@
+#ifndef STORM_STORAGE_BISIMULATION_DETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+#define STORM_STORAGE_BISIMULATION_DETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+
+#include "src/storage/bisimulation/BisimulationDecomposition.h"
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+namespace storm {
+    namespace utility {
+        template <typename ValueType> class ConstantsComparator;
+    }
+    
+    namespace storage {
+        
+        /*!
+         * This class represents the decomposition of a deterministic model into its bisimulation quotient.
+         */
+        template<typename ModelType>
+        class DeterministicModelBisimulationDecomposition : public BisimulationDecomposition<ModelType, bisimulation::DeterministicBlockData> {
+        public:
+            typedef bisimulation::DeterministicBlockData BlockDataType;
+            typedef typename ModelType::ValueType ValueType;
+            typedef typename ModelType::RewardModelType RewardModelType;
+            
+            /*!
+             * Computes the bisimulation relation for the given model. Which kind of bisimulation is computed, is
+             * customizable via the given options.
+             *
+             * @param model The model to decompose.
+             * @param options The options that customize the computed bisimulation.
+             */
+            DeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, BlockDataType>::Options const& options = typename BisimulationDecomposition<ModelType, BlockDataType>::Options());
+            
+        protected:
+            virtual std::pair<storm::storage::BitVector, storm::storage::BitVector> getStatesWithProbability01() override;
+
+            virtual void initializeMeasureDrivenPartition() override;
+            
+            virtual void initializeLabelBasedPartition() override;
+            
+            virtual void buildQuotient() override;
+            
+            virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) override;
+
+        private:
+            // Refines the predecessor blocks wrt. strong bisimulation.
+            void refinePredecessorBlocksOfSplitterStrong(std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+
+            /*!
+             * Performs the necessary steps to compute a weak bisimulation on a DTMC.
+             */
+            void initializeWeakDtmcBisimulation();
+            
+            /*!
+             * Splits all blocks of the current partition into a block that contains all divergent states and another
+             * block containing the non-divergent states.
+             */
+            void splitOffDivergentStates();
+            
+            /*!
+             * Initializes the vector of silent probabilities.
+             */
+            void initializeSilentProbabilities();
+            
+            // Retrieves the probability of going into the splitter for the given state.
+            ValueType const& getProbabilityToSplitter(storm::storage::sparse::state_type const& state) const;
+            
+            // Retrieves the silent probability for the given state.
+            ValueType getSilentProbability(storm::storage::sparse::state_type const& state) const;
+            
+            // Retrieves whether the given state is silent.
+            bool isSilent(storm::storage::sparse::state_type const& state) const;
+            
+            // Retrieves whether the given state has a non-zero silent probability.
+            bool hasNonZeroSilentProbability(storm::storage::sparse::state_type const& state) const;
+            
+            // Retrieves whether the given predecessor of the splitters possibly needs refinement.
+            bool possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& predecessorBlock) const;
+            
+            // Moves the given state to the position marked by marker1 moves the marker one step further.
+            void moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock);
+
+            // Moves the given state to the position marked by marker2 the marker one step further.
+            void moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock);
+            
+            // Moves the given state to a proper place in the splitter, depending on where the predecessor is located.
+            void moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip);
+            
+            // Increases the probability of moving to the current splitter for the given state.
+            void increaseProbabilityToSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType> const& predecessorBlock, ValueType const& value);
+            
+            // Explores the remaining states of the splitter.
+            void exploreRemainingStatesOfSplitter(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks);
+            
+            // Updates the silent probabilities of the states in the block based on the probabilities of going to the splitter.
+            void updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(bisimulation::Block<BlockDataType>& block);
+
+            // Updates the silent probabilities of the states in the block based on a forward exploration of the transitions
+            // of the states.
+            void updateSilentProbabilitiesBasedOnTransitions(bisimulation::Block<BlockDataType>& block);
+            
+            // Refines the predecessor blocks of the splitter wrt. weak bisimulation in DTMCs.
+            void refinePredecessorBlocksOfSplitterWeak(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            // Refines the given block wrt to weak bisimulation in DTMCs.
+            void refinePredecessorBlockOfSplitterWeak(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            // Converts the one-step probabilities of going into the splitter into the conditional probabilities needed
+            // for weak bisimulation (on DTMCs).
+            void computeConditionalProbabilitiesForNonSilentStates(bisimulation::Block<BlockDataType>& block);
+            
+            // Computes the (indices of the) blocks of non-silent states within the block.
+            std::vector<uint_fast64_t> computeNonSilentBlocks(bisimulation::Block<BlockDataType>& block);
+            
+            // Computes a labeling for all states of the block that identifies in which block they need to end up.
+            std::vector<storm::storage::BitVector> computeWeakStateLabelingBasedOnNonSilentBlocks(bisimulation::Block<BlockDataType> const& block, std::vector<uint_fast64_t> const& nonSilentBlockIndices);
+            
+            // Inserts the block into the list of predecessors if it is not already contained.
+            void insertIntoPredecessorList(bisimulation::Block<BlockDataType>& predecessorBlock, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks);
+            
+            // A vector that holds the probabilities of states going into the splitter. This is used by the method that
+            // refines a block based on probabilities.
+            std::vector<ValueType> probabilitiesToCurrentSplitter;
+            
+            // A vector mapping each state to its silent probability.
+            std::vector<ValueType> silentProbabilities;
+        };
+    }
+}
+
+#endif /* STORM_STORAGE_BISIMULATION_DETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ */
\ No newline at end of file
diff --git a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp
new file mode 100644
index 000000000..21dd36d89
--- /dev/null
+++ b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp
@@ -0,0 +1,425 @@
+#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
+
+#include "src/models/sparse/Mdp.h"
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/utility/graph.h"
+
+#include "src/utility/macros.h"
+#include "src/exceptions/IllegalFunctionCallException.h"
+
+#include "src/adapters/CarlAdapter.h"
+
+namespace storm {
+    namespace storage {
+        
+        using namespace bisimulation;
+        
+        template<typename ModelType>
+        NondeterministicModelBisimulationDecomposition<ModelType>::NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>::Options const& options) : BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>(model, model.getTransitionMatrix().transpose(false), options), choiceToStateMapping(model.getNumberOfChoices()), quotientDistributions(model.getNumberOfChoices()), orderedQuotientDistributions(model.getNumberOfChoices()) {
+            STORM_LOG_THROW(options.type == BisimulationType::Strong, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation is currently not supported for nondeterministic models.");
+        }
+        
+        template<typename ModelType>
+        std::pair<storm::storage::BitVector, storm::storage::BitVector> NondeterministicModelBisimulationDecomposition<ModelType>::getStatesWithProbability01() {
+            STORM_LOG_THROW(static_cast<bool>(this->options.optimalityType), storm::exceptions::IllegalFunctionCallException, "Can only compute states with probability 0/1 with an optimization direction (min/max).");
+            if (this->options.optimalityType.get() == OptimizationDirection::Minimize) {
+                return storm::utility::graph::performProb01Min(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(), this->model.getBackwardTransitions(), this->options.phiStates.get(), this->options.psiStates.get());
+            } else {
+                return storm::utility::graph::performProb01Max(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(), this->model.getBackwardTransitions(), this->options.phiStates.get(), this->options.psiStates.get());
+            }
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::initialize() {
+            this->createChoiceToStateMapping();
+            this->initializeQuotientDistributions();
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::createChoiceToStateMapping() {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (storm::storage::sparse::state_type state = 0; state < this->model.getNumberOfStates(); ++state) {
+                for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                    choiceToStateMapping[choice] = state;
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::initializeQuotientDistributions() {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            
+            for (auto const& block : this->partition.getBlocks()) {
+                if (block->data().absorbing()) {
+                    // If the block is marked as absorbing, we need to create the corresponding distributions.
+                    for (auto stateIt = this->partition.begin(*block), stateIte = this->partition.end(*block); stateIt != stateIte; ++stateIt) {
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[*stateIt]; choice < nondeterministicChoiceIndices[*stateIt + 1]; ++choice) {
+                            this->quotientDistributions[choice].addProbability(block->getId(), storm::utility::one<ValueType>());
+                            orderedQuotientDistributions[choice] = &this->quotientDistributions[choice];
+                        }
+                    }
+                } else {
+                    // Otherwise, we compute the probabilities from the transition matrix.
+                    for (auto stateIt = this->partition.begin(*block), stateIte = this->partition.end(*block); stateIt != stateIte; ++stateIt) {
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[*stateIt]; choice < nondeterministicChoiceIndices[*stateIt + 1]; ++choice) {
+                            for (auto entry : this->model.getTransitionMatrix().getRow(choice)) {
+                                if (!this->comparator.isZero(entry.getValue())) {
+                                    this->quotientDistributions[choice].addProbability(this->partition.getBlock(entry.getColumn()).getId(), entry.getValue());
+                                }
+                            }
+                            orderedQuotientDistributions[choice] = &this->quotientDistributions[choice];
+                        }
+                    }
+                }
+            }
+            
+            for (auto state = 0; state < this->model.getNumberOfStates(); ++state) {
+                updateOrderedQuotientDistributions(state);
+            }
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::updateOrderedQuotientDistributions(storm::storage::sparse::state_type state) {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            std::sort(this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state], this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state + 1],
+                      [this] (storm::storage::Distribution<ValueType> const* dist1, storm::storage::Distribution<ValueType> const* dist2) {
+                          return dist1->less(*dist2, this->comparator);
+                      });
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
+            // In order to create the quotient model, we need to construct
+            // (a) the new transition matrix,
+            // (b) the new labeling,
+            // (c) the new reward structures.
+            
+            // Prepare a matrix builder for (a).
+            storm::storage::SparseMatrixBuilder<ValueType> builder(0, this->size(), 0, false, true, this->size());
+            
+            // Prepare the new state labeling for (b).
+            storm::models::sparse::StateLabeling newLabeling(this->size());
+            std::set<std::string> atomicPropositionsSet = this->options.respectedAtomicPropositions.get();
+            atomicPropositionsSet.insert("init");
+            std::vector<std::string> atomicPropositions = std::vector<std::string>(atomicPropositionsSet.begin(), atomicPropositionsSet.end());
+            for (auto const& ap : atomicPropositions) {
+                newLabeling.addLabel(ap);
+            }
+            
+            // If the model had state rewards, we need to build the state rewards for the quotient as well.
+            boost::optional<std::vector<ValueType>> stateRewards;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                stateRewards = std::vector<ValueType>(this->blocks.size());
+            }
+            
+            // Now build (a) and (b) by traversing all blocks.
+            uint_fast64_t currentRow = 0;
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) {
+                auto const& block = this->blocks[blockIndex];
+                
+                // Open new row group for the new meta state.
+                builder.newRowGroup(currentRow);
+                
+                // Pick one representative state. For strong bisimulation it doesn't matter which state it is, because
+                // they all behave equally.
+                storm::storage::sparse::state_type representativeState = *block.begin();
+                Block<BlockDataType> const& oldBlock = this->partition.getBlock(representativeState);
+                
+                // If the block is absorbing, we simply add a self-loop.
+                if (oldBlock.data().absorbing()) {
+                    builder.addNextValue(currentRow, blockIndex, storm::utility::one<ValueType>());
+                    ++currentRow;
+                    
+                    // If the block has a special representative state, we retrieve it now.
+                    if (oldBlock.data().hasRepresentativeState()) {
+                        representativeState = oldBlock.data().representativeState();
+                    }
+                    
+                    // Add all of the selected atomic propositions that hold in the representative state to the state
+                    // representing the block.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                } else {
+                    // Add the outgoing choices of the block.
+                    for (uint_fast64_t choice = nondeterministicChoiceIndices[representativeState]; choice < nondeterministicChoiceIndices[representativeState + 1]; ++choice) {
+                        // If the choice is the same as the last one, we do not need to add it.
+                        if (choice > nondeterministicChoiceIndices[representativeState] && quotientDistributions[choice - 1].equals(quotientDistributions[choice], this->comparator)) {
+                            continue;
+                        }
+                        
+                        for (auto entry : quotientDistributions[choice]) {
+                            builder.addNextValue(currentRow, entry.first, entry.second);
+                        }
+                        ++currentRow;
+                    }
+                    
+                    // Otherwise add all atomic propositions to the equivalence class that the representative state
+                    // satisfies.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                }
+                
+                // If the model has state rewards, we simply copy the state reward of the representative state, because
+                // all states in a block are guaranteed to have the same state reward.
+                if (this->options.keepRewards && this->model.hasRewardModel()) {
+                    typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                    stateRewards.get()[blockIndex] = nameRewardModelPair->second.getStateRewardVector()[representativeState];
+                }
+            }
+            
+            // Now check which of the blocks of the partition contain at least one initial state.
+            for (auto initialState : this->model.getInitialStates()) {
+                Block<BlockDataType> const& initialBlock = this->partition.getBlock(initialState);
+                newLabeling.addLabelToState("init", initialBlock.getId());
+            }
+            
+            // Construct the reward model mapping.
+            std::unordered_map<std::string, typename ModelType::RewardModelType> rewardModels;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                rewardModels.insert(std::make_pair(nameRewardModelPair->first, typename ModelType::RewardModelType(stateRewards)));
+            }
+            
+            // Finally construct the quotient model.
+            this->quotient = std::shared_ptr<ModelType>(new ModelType(builder.build(), std::move(newLabeling), std::move(rewardModels)));
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& block) const {
+            return block.getNumberOfStates() > 1 && !block.data().absorbing();
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::updateQuotientDistributionsOfPredecessors(Block<BlockDataType> const& newBlock, Block<BlockDataType> const& oldBlock, std::deque<Block<BlockDataType>*>& splitterQueue) {
+            uint_fast64_t lastState = 0;
+            bool lastStateInitialized = false;
+            
+            for (auto stateIt = this->partition.begin(newBlock), stateIte = this->partition.end(newBlock); stateIt != stateIte; ++stateIt) {
+                for (auto predecessorEntry : this->backwardTransitions.getRow(*stateIt)) {
+                    if (this->comparator.isZero(predecessorEntry.getValue())) {
+                        continue;
+                    }
+                    
+                    storm::storage::sparse::state_type predecessorChoice = predecessorEntry.getColumn();
+                    storm::storage::sparse::state_type predecessorState = choiceToStateMapping[predecessorChoice];
+                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessorState);
+                    
+                    // If the predecessor block is marked as absorbing, we do not need to update anything.
+                    if (predecessorBlock.data().absorbing()) {
+                        continue;
+                    }
+                    
+                    // If the predecessor block is not marked as to-refined, we do so now.
+                    if (!predecessorBlock.data().splitter()) {
+                        predecessorBlock.data().setSplitter();
+                        splitterQueue.push_back(&predecessorBlock);
+                    }
+                    
+                    if (lastStateInitialized) {
+                        // If we have skipped to the choices of the next state, we need to repair the order of the
+                        // distributions for the last state.
+                        if (lastState != predecessorState) {
+                            updateOrderedQuotientDistributions(lastState);
+                            lastState = predecessorState;
+                        }
+                    } else {
+                        lastStateInitialized = true;
+                        lastState = choiceToStateMapping[predecessorChoice];
+                    }
+                    
+                    // Now shift the probability from this transition from the old block to the new one.
+                    this->quotientDistributions[predecessorChoice].shiftProbability(oldBlock.getId(), newBlock.getId(), predecessorEntry.getValue());
+                }
+            }
+            
+            if (lastStateInitialized) {
+                updateOrderedQuotientDistributions(lastState);
+            }
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkQuotientDistributions() const {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (auto state = 0; state < this->model.getNumberOfStates(); ++state) {
+                for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                    storm::storage::Distribution<ValueType> distribution;
+                    for (auto const& element : this->model.getTransitionMatrix().getRow(choice)) {
+                        distribution.addProbability(this->partition.getBlock(element.getColumn()).getId(), element.getValue());
+                    }
+
+                    if (!distribution.equals(quotientDistributions[choice])) {
+                        std::cout << "the distributions for choice " << choice << " of state " << state << " do not match." << std::endl;
+                        std::cout << "is: " << quotientDistributions[choice] << " but should be " << distribution << std::endl;
+                        exit(-1);
+                    }
+
+                    bool less1 = distribution.less(quotientDistributions[choice], this->comparator);
+                    bool less2 = quotientDistributions[choice].less(distribution, this->comparator);
+                    
+                    if (distribution.equals(quotientDistributions[choice]) && (less1 || less2)) {
+                        std::cout << "mismatch of equality and less for " << std::endl;
+                        std::cout << quotientDistributions[choice] << " vs " << distribution << std::endl;
+                        exit(-1);
+                    }
+                }
+                
+                for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1] - 1; ++choice) {
+                    if (orderedQuotientDistributions[choice + 1]->less(*orderedQuotientDistributions[choice], this->comparator)) {
+                        std::cout << "choice " << (choice+1) << " is less than predecessor" << std::endl;
+                        std::cout << *orderedQuotientDistributions[choice] << " should be less than " << *orderedQuotientDistributions[choice + 1] << std::endl;
+                        exit(-1);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::printDistributions(uint_fast64_t state) const {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                std::cout << quotientDistributions[choice] << std::endl;
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkBlockStable(bisimulation::Block<BlockDataType> const& newBlock) const {
+            std::cout << "checking stability of new block " << newBlock.getId() << " of size " << newBlock.getNumberOfStates() << std::endl;
+            for (auto stateIt1 = this->partition.begin(newBlock), stateIte1 = this->partition.end(newBlock); stateIt1 != stateIte1; ++stateIt1) {
+                for (auto stateIt2 = this->partition.begin(newBlock), stateIte2 = this->partition.end(newBlock); stateIt2 != stateIte2; ++stateIt2) {
+                    bool less1 = quotientDistributionsLess(*stateIt1, *stateIt2);
+                    bool less2 = quotientDistributionsLess(*stateIt2, *stateIt1);
+                    if (less1 || less2) {
+                        std::cout << "the partition is not stable for the states " << *stateIt1 << " and " << *stateIt2 << std::endl;
+                        std::cout << "less1 " << less1 << " and less2 " << less2 << std::endl;
+                        
+                        std::cout << "distributions of state " << *stateIt1 << std::endl;
+                        this->printDistributions(*stateIt1);
+                        std::cout << "distributions of state " << *stateIt2 << std::endl;
+                        this->printDistributions(*stateIt2);
+                        exit(-1);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkDistributionsDifferent(bisimulation::Block<BlockDataType> const& block, storm::storage::sparse::state_type end) const {
+            for (auto stateIt1 = this->partition.begin(block), stateIte1 = this->partition.end(block); stateIt1 != stateIte1; ++stateIt1) {
+                for (auto stateIt2 = this->partition.begin() + block.getEndIndex(), stateIte2 = this->partition.begin() + end; stateIt2 != stateIte2; ++stateIt2) {
+                    if (!quotientDistributionsLess(*stateIt1, *stateIt2)) {
+                        std::cout << "distributions are not less, even though they should be!" << std::endl;
+                        exit(-3);
+                    } else {
+                        std::cout << "less:" << std::endl;
+                        this->printDistributions(*stateIt1);
+                        std::cout << "and" << std::endl;
+                        this->printDistributions(*stateIt2);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::splitBlockAccordingToCurrentQuotientDistributions(Block<BlockDataType>& block, std::deque<Block<BlockDataType>*>& splitterQueue) {
+            std::list<Block<BlockDataType>*> newBlocks;
+            bool split = this->partition.splitBlock(block,
+                                                    [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                                                        bool result = quotientDistributionsLess(state1, state2);
+//                                                        std::cout << state1 << " is " << (!result ? "not" : "") << " smaller than " << state2 << std::endl;
+                                                        return result;
+                                                    },
+                                                    [this, &block, &splitterQueue, &newBlocks] (Block<BlockDataType>& newBlock) {
+                                                        newBlocks.push_back(&newBlock);
+                                                        
+//                                                        this->checkBlockStable(newBlock);
+//                                                        this->checkDistributionsDifferent(block, block.getEndIndex());
+//                                                        this->checkQuotientDistributions();
+                                                    });
+            
+            // The quotient distributions of the predecessors of block do not need to be updated, since the probability
+            // will go to the block with the same id as before.
+            
+//            std::cout << "partition after split: " << std::endl;
+//            this->partition.print();
+            
+            // defer updating the quotient distributions until *after* all splits, because
+            // it otherwise influences the subsequent splits!
+            for (auto el : newBlocks) {
+                this->updateQuotientDistributionsOfPredecessors(*el, block, splitterQueue);
+            }
+            
+//            this->checkQuotientDistributions();
+            
+            return split;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) const {
+            STORM_LOG_TRACE("Comparing the quotient distributions of state " << state1 << " and " << state2 << ".");
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+
+//            std::cout << "distributions of state " << state1 << std::endl;
+//            this->printDistributions(state1);
+//            std::cout << "distributions of state " << state2 << std::endl;
+//            this->printDistributions(state2);
+            
+            auto firstIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1];
+            auto firstIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1 + 1];
+            auto secondIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2];
+            auto secondIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2 + 1];
+            
+            for (; firstIt != firstIte && secondIt != secondIte; ++firstIt, ++secondIt) {
+                // If the current distributions are in a less-than relationship, we can return a result.
+                if ((*firstIt)->less(**secondIt, this->comparator)) {
+                    return true;
+                } else if ((*secondIt)->less(**firstIt, this->comparator)) {
+                    return false;
+                }
+                
+                // If the distributions matched, we need to advance both distribution iterators to the next distribution
+                // that is larger.
+                while (firstIt != firstIte && std::next(firstIt) != firstIte && !(*firstIt)->less(**std::next(firstIt), this->comparator)) {
+                    ++firstIt;
+                }
+                while (secondIt != secondIte && std::next(secondIt) != secondIte && !(*secondIt)->less(**std::next(secondIt), this->comparator)) {
+                    ++secondIt;
+                }
+            }
+            
+            if (firstIt == firstIte && secondIt != secondIte) {
+                return true;
+            }
+            return false;
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            if (!possiblyNeedsRefinement(splitter)) {
+                return;
+            }
+            
+            STORM_LOG_TRACE("Refining block " << splitter.getId());
+            
+            splitBlockAccordingToCurrentQuotientDistributions(splitter, splitterQueue);
+        }
+        
+        template class NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>;
+        
+#ifdef STORM_HAVE_CARL
+        template class NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<storm::RationalFunction>>;
+#endif
+        
+    }
+}
diff --git a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h
new file mode 100644
index 000000000..18a6e8e58
--- /dev/null
+++ b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h
@@ -0,0 +1,84 @@
+#ifndef STORM_STORAGE_BISIMULATION_NONDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+#define STORM_STORAGE_BISIMULATION_NONDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+
+#include "src/storage/bisimulation/BisimulationDecomposition.h"
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+#include "src/storage/Distribution.h"
+
+namespace storm {
+    namespace utility {
+        template <typename ValueType> class ConstantsComparator;
+    }
+    
+    namespace storage {
+        
+        /*!
+         * This class represents the decomposition of a nondeterministic model into its bisimulation quotient.
+         */
+        template<typename ModelType>
+        class NondeterministicModelBisimulationDecomposition : public BisimulationDecomposition<ModelType, bisimulation::DeterministicBlockData> {
+        public:
+            typedef bisimulation::DeterministicBlockData BlockDataType;
+            typedef typename ModelType::ValueType ValueType;
+            typedef typename ModelType::RewardModelType RewardModelType;
+            
+            /*!
+             * Computes the bisimulation relation for the given model. Which kind of bisimulation is computed, is
+             * customizable via the given options.
+             *
+             * @param model The model to decompose.
+             * @param options The options that customize the computed bisimulation.
+             */
+            NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, BlockDataType>::Options const& options = typename BisimulationDecomposition<ModelType, BlockDataType>::Options());
+            
+        protected:
+            virtual std::pair<storm::storage::BitVector, storm::storage::BitVector> getStatesWithProbability01() override;
+            
+            virtual void buildQuotient() override;
+            
+            virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) override;
+            
+            virtual void initialize() override;
+            
+        private:
+            // Creates the mapping from the choice indices to the states.
+            void createChoiceToStateMapping();
+            
+            // Initializes the quotient distributions wrt. to the current partition.
+            void initializeQuotientDistributions();
+            
+            // Retrieves whether the given block possibly needs refinement.
+            bool possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& block) const;
+            
+            // Splits the given block according to the current quotient distributions.
+            bool splitBlockAccordingToCurrentQuotientDistributions(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            // Retrieves whether the quotient distributions of state 1 are considered to be less than the ones of state 2.
+            bool quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) const;
+            
+            // Updates the ordered list of quotient distribution for the given state.
+            void updateOrderedQuotientDistributions(storm::storage::sparse::state_type state);
+            
+            // Updates the quotient distributions of the predecessors of the new block by taking the probability mass
+            // away from the old block.
+            void updateQuotientDistributionsOfPredecessors(bisimulation::Block<BlockDataType> const& newBlock, bisimulation::Block<BlockDataType> const& oldBlock, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            bool checkQuotientDistributions() const;
+            bool checkBlockStable(bisimulation::Block<BlockDataType> const& newBlock) const;
+            bool printDistributions(storm::storage::sparse::state_type state) const;
+            bool checkDistributionsDifferent(bisimulation::Block<BlockDataType> const& block, storm::storage::sparse::state_type end) const;
+            
+            // A mapping from choice indices to the state state that has this choice.
+            std::vector<storm::storage::sparse::state_type> choiceToStateMapping;
+            
+            // A vector that holds the quotient distributions for all nondeterministic choices of all states.
+            std::vector<storm::storage::Distribution<ValueType>> quotientDistributions;
+            
+            // A vector that stores for each state the ordered list of quotient distributions.
+            std::vector<storm::storage::Distribution<ValueType> const*> orderedQuotientDistributions;
+        };
+    }
+}
+
+#endif /* STORM_STORAGE_BISIMULATION_NONDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ */
\ No newline at end of file
diff --git a/src/storage/bisimulation/Partition.cpp b/src/storage/bisimulation/Partition.cpp
new file mode 100644
index 000000000..698eb5cd7
--- /dev/null
+++ b/src/storage/bisimulation/Partition.cpp
@@ -0,0 +1,379 @@
+#include "src/storage/bisimulation/Partition.h"
+
+#include <iostream>
+
+#include "src/storage/bisimulation/DeterministicBlockData.h"
+
+#include "src/utility/macros.h"
+#include "src/exceptions/InvalidArgumentException.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+            template<typename DataType>
+            Partition<DataType>::Partition(std::size_t numberOfStates) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates) {
+                blocks.emplace_back(new Block<DataType>(0, numberOfStates, nullptr, nullptr, blocks.size()));
+                
+                // Set up the different parts of the internal structure.
+                for (storm::storage::sparse::state_type state = 0; state < numberOfStates; ++state) {
+                    states[state] = state;
+                    positions[state] = state;
+                    stateToBlockMapping[state] = blocks.back().get();
+                }
+            }
+            
+            template<typename DataType>
+            Partition<DataType>::Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, boost::optional<storm::storage::sparse::state_type> representativeProb1State) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates) {
+                storm::storage::sparse::state_type position = 0;
+                Block<DataType>* firstBlock = nullptr;
+                Block<DataType>* secondBlock = nullptr;
+                Block<DataType>* thirdBlock = nullptr;
+                if (!prob0States.empty()) {
+                    blocks.emplace_back(new Block<DataType>(0, prob0States.getNumberOfSetBits(), nullptr, nullptr, blocks.size()));
+                    firstBlock = blocks.front().get();
+                    
+                    for (auto state : prob0States) {
+                        states[position] = state;
+                        positions[state] = position;
+                        stateToBlockMapping[state] = firstBlock;
+                        ++position;
+                    }
+                    firstBlock->data().setAbsorbing(true);
+                }
+                
+                if (!prob1States.empty()) {
+                    blocks.emplace_back(new Block<DataType>(position, position + prob1States.getNumberOfSetBits(), firstBlock, nullptr, blocks.size()));
+                    secondBlock = blocks.back().get();
+                    
+                    for (auto state : prob1States) {
+                        states[position] = state;
+                        positions[state] = position;
+                        stateToBlockMapping[state] = secondBlock;
+                        ++position;
+                    }
+                    secondBlock->data().setAbsorbing(true);
+                    secondBlock->data().setRepresentativeState(representativeProb1State.get());
+                }
+                
+                storm::storage::BitVector otherStates = ~(prob0States | prob1States);
+                if (!otherStates.empty()) {
+                    blocks.emplace_back(new Block<DataType>(position, numberOfStates, secondBlock, nullptr, blocks.size()));
+                    thirdBlock = blocks.back().get();
+                    
+                    for (auto state : otherStates) {
+                        states[position] = state;
+                        positions[state] = position;
+                        stateToBlockMapping[state] = thirdBlock;
+                        ++position;
+                    }
+                }
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                std::swap(this->states[this->positions[state1]], this->states[this->positions[state2]]);
+                std::swap(this->positions[state1], this->positions[state2]);
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2) {
+                storm::storage::sparse::state_type state1 = this->states[position1];
+                storm::storage::sparse::state_type state2 = this->states[position2];
+                std::swap(this->states[position1], this->states[position2]);
+                this->positions[state1] = position2;
+                this->positions[state2] = position1;
+            }
+            
+            template<typename DataType>
+            storm::storage::sparse::state_type const& Partition<DataType>::getPosition(storm::storage::sparse::state_type state) const {
+                return this->positions[state];
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position) {
+                this->positions[state] = position;
+            }
+            
+            template<typename DataType>
+            storm::storage::sparse::state_type const& Partition<DataType>::getState(storm::storage::sparse::state_type position) const {
+                return this->states[position];
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::mapStatesToBlock(Block<DataType>& block, std::vector<storm::storage::sparse::state_type>::iterator first, std::vector<storm::storage::sparse::state_type>::iterator last) {
+                for (; first != last; ++first) {
+                    this->stateToBlockMapping[*first] = &block;
+                }
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::mapStatesToPositions(std::vector<storm::storage::sparse::state_type>::const_iterator first, std::vector<storm::storage::sparse::state_type>::const_iterator last) {
+                storm::storage::sparse::state_type position = std::distance(this->states.cbegin(), first);
+                for (; first != last; ++first, ++position) {
+                    this->positions[*first] = position;
+                }
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::mapStatesToPositions(Block<DataType> const& block) {
+                mapStatesToPositions(this->begin(block), this->end(block));
+            }
+            
+            template<typename DataType>
+            Block<DataType>& Partition<DataType>::getBlock(storm::storage::sparse::state_type state) {
+                return *this->stateToBlockMapping[state];
+            }
+            
+            template<typename DataType>
+            Block<DataType> const& Partition<DataType>::getBlock(storm::storage::sparse::state_type state) const {
+                return *this->stateToBlockMapping[state];
+            }
+
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::iterator Partition<DataType>::begin(Block<DataType> const& block) {
+                auto it = this->states.begin();
+                std::advance(it, block.getBeginIndex());
+                return it;
+            }
+            
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::const_iterator Partition<DataType>::begin(Block<DataType> const& block) const {
+                auto it = this->states.begin();
+                std::advance(it, block.getBeginIndex());
+                return it;
+            }
+
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::iterator Partition<DataType>::end(Block<DataType> const& block) {
+                auto it = this->states.begin();
+                std::advance(it, block.getEndIndex());
+                return it;
+            }
+
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::const_iterator Partition<DataType>::end(Block<DataType> const& block) const {
+                auto it = this->states.begin();
+                std::advance(it, block.getEndIndex());
+                return it;
+            }
+            
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::iterator Partition<DataType>::begin() {
+                return this->states.begin();
+            }
+            
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::const_iterator Partition<DataType>::begin() const {
+                return this->states.begin();
+            }
+            
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::iterator Partition<DataType>::end() {
+                return this->states.end();
+            }
+            
+            template<typename DataType>
+            std::vector<storm::storage::sparse::state_type>::const_iterator Partition<DataType>::end() const {
+                return this->states.end();
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::sortRange(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions) {
+                std::sort(this->states.begin() + beginIndex, this->states.begin() + endIndex, less);
+                
+                if (updatePositions) {
+                    mapStatesToPositions(this->states.begin() + beginIndex, this->states.begin() + endIndex);
+                }
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::sortBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions) {
+                sortRange(block.getBeginIndex(), block.getEndIndex(), less, updatePositions);
+            }
+            
+            template<typename DataType>
+            std::vector<uint_fast64_t> Partition<DataType>::computeRangesOfEqualValue(uint_fast64_t startIndex, uint_fast64_t endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less) {
+                auto it = this->states.cbegin() + startIndex;
+                auto ite = this->states.cbegin() + endIndex;
+                
+                std::vector<storm::storage::sparse::state_type>::const_iterator upperBound;
+                std::vector<uint_fast64_t> result;
+                result.push_back(startIndex);
+                do {
+                    upperBound = std::upper_bound(it, ite, *it, less);
+                    result.push_back(std::distance(this->states.cbegin(), upperBound));
+                    it = upperBound;
+                } while (upperBound != ite);
+                
+                return result;
+            }
+            
+            template<typename DataType>
+            std::pair<typename std::vector<std::unique_ptr<Block<DataType>>>::iterator, bool> Partition<DataType>::splitBlock(Block<DataType>& block, storm::storage::sparse::state_type position) {
+                STORM_LOG_THROW(position >= block.getBeginIndex() && position <= block.getEndIndex(), storm::exceptions::InvalidArgumentException, "Cannot split block at illegal position.");
+                STORM_LOG_TRACE("Splitting " << block.getId() << " at position " << position << " (begin was " << block.getBeginIndex() << ".");
+
+                // In case one of the resulting blocks would be empty, we simply return the current block and do not create
+                // a new one.
+                if (position == block.getBeginIndex() || position == block.getEndIndex()) {
+                    auto it = blocks.begin();
+                    std::advance(it, block.getId());
+                    return std::make_pair(it, false);
+                }
+                
+                // Actually create the new block.
+                blocks.emplace_back(new Block<DataType>(block.getBeginIndex(), position, block.getPreviousBlockPointer(), &block, blocks.size()));
+                auto newBlockIt = std::prev(blocks.end());
+                
+                // Resize the current block appropriately.
+                block.setBeginIndex(position);
+                
+                // Update the mapping of the states in the newly created block.
+                this->mapStatesToBlock(**newBlockIt, this->begin(**newBlockIt), this->end(**newBlockIt));
+                
+                return std::make_pair(newBlockIt, true);
+            }
+            
+            template<typename DataType>
+            bool Partition<DataType>::splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback) {
+                // Sort the block, but leave the positions untouched.
+                this->sortBlock(block, less, false);
+                
+                auto originalBegin = block.getBeginIndex();
+                auto originalEnd = block.getEndIndex();
+                
+                auto it = this->states.cbegin() + block.getBeginIndex();
+                auto ite = this->states.cbegin() + block.getEndIndex();
+                
+                bool wasSplit = false;
+                std::vector<storm::storage::sparse::state_type>::const_iterator upperBound;
+                do {
+                    upperBound = std::upper_bound(it, ite, *it, less);
+                    
+                    if (upperBound != ite) {
+                        wasSplit = true;
+                        auto result = this->splitBlock(block, std::distance(this->states.cbegin(), upperBound));
+                        newBlockCallback(**result.first);
+                    }
+                    it = upperBound;
+                } while (upperBound != ite);
+                
+                // Finally, repair the positions mapping.
+                mapStatesToPositions(this->states.begin() + originalBegin, this->states.begin() + originalEnd);
+
+                return wasSplit;
+            }
+            
+            // Splits the block by sorting the states according to the given function and then identifying the split
+            // points.
+            template<typename DataType>
+            bool Partition<DataType>::splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less) {
+                return this->splitBlock(block, less, [] (Block<DataType>& block) {});
+            }
+            
+            template<typename DataType>
+            bool Partition<DataType>::split(std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback) {
+                bool result = false;
+                // Since the underlying storage of the blocks may change during iteration, we remember the current size
+                // and iterate over indices. This assumes that new blocks will be added at the end of the blocks vector.
+                std::size_t currentSize = this->size();
+                for (uint_fast64_t index = 0; index < currentSize; ++index) {
+                    result |= splitBlock(*blocks[index], less, newBlockCallback);
+                }
+                return result;
+            }
+            
+            template<typename DataType>
+            bool Partition<DataType>::split(std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less) {
+                return this->split(less, [] (Block<DataType>& block) {});
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::splitStates(Block<DataType>& block, storm::storage::BitVector const& states) {
+                this->splitBlock(block, [&states] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return states.get(a) && !states.get(b); });
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::splitStates(storm::storage::BitVector const& states) {
+                this->split([&states] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return states.get(a) && !states.get(b); });
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::sortBlock(Block<DataType> const& block) {
+                std::sort(this->begin(block), this->end(block), [] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return a < b; });
+                mapStatesToPositions(block);
+            }
+            
+            template<typename DataType>
+            Block<DataType>& Partition<DataType>::insertBlock(Block<DataType>& block) {
+                // Find the beginning of the new block.
+                storm::storage::sparse::state_type begin = block.hasPreviousBlock() ? block.getPreviousBlock().getEndIndex() : 0;
+                
+                // Actually insert the new block.
+                blocks.emplace_back(new Block<DataType>(begin, block.getBeginIndex(), block.getPreviousBlockPointer(), &block, blocks.size()));
+                Block<DataType>& newBlock = *blocks.back();
+                
+                // Update the mapping of the states in the newly created block.
+                for (auto it = this->begin(newBlock), ite = this->end(newBlock); it != ite; ++it) {
+                    stateToBlockMapping[*it] = &newBlock;
+                }
+                
+                return newBlock;
+            }
+            
+            template<typename DataType>
+            std::vector<std::unique_ptr<Block<DataType>>> const& Partition<DataType>::getBlocks() const {
+                return this->blocks;
+            }
+            
+            template<typename DataType>
+            std::vector<std::unique_ptr<Block<DataType>>>& Partition<DataType>::getBlocks() {
+                return this->blocks;
+            }
+            
+            template<typename DataType>
+            bool Partition<DataType>::check() const {
+                for (uint_fast64_t state = 0; state < this->positions.size(); ++state) {
+                    STORM_LOG_ASSERT(this->states[this->positions[state]] == state, "Position mapping corrupted.");
+                }
+                for (auto const& blockPtr : this->blocks) {
+                    STORM_LOG_ASSERT(blockPtr->check(), "Block corrupted.");
+                    for (auto stateIt = this->begin(*blockPtr), stateIte = this->end(*blockPtr); stateIt != stateIte; ++stateIt) {
+                        STORM_LOG_ASSERT(this->stateToBlockMapping[*stateIt] == blockPtr.get(), "Block mapping corrupted.");
+                    }
+                }
+                return true;
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::print() const {
+                for (auto const& block : this->blocks) {
+                    block->print(*this);
+                }
+                std::cout << "states in partition" << std::endl;
+                for (auto const& state : states) {
+                    std::cout << state << " ";
+                }
+                std::cout << std::endl << "positions: " << std::endl;
+                for (auto const& index : positions) {
+                    std::cout << index << " ";
+                }
+                std::cout << std::endl << "state to block mapping: " << std::endl;
+                for (auto const& block : stateToBlockMapping) {
+                    std::cout << block << "[" << block->getId() <<"] ";
+                }
+                std::cout << std::endl;
+                std::cout << "size: " << blocks.size() << std::endl;
+                STORM_LOG_ASSERT(this->check(), "Partition corrupted.");
+            }
+            
+            template<typename DataType>
+            std::size_t Partition<DataType>::size() const {
+                return blocks.size();
+            }
+
+            template class Partition<DeterministicBlockData>;
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/storage/bisimulation/Partition.h b/src/storage/bisimulation/Partition.h
new file mode 100644
index 000000000..b19830a44
--- /dev/null
+++ b/src/storage/bisimulation/Partition.h
@@ -0,0 +1,181 @@
+#ifndef STORM_STORAGE_BISIMULATION_PARTITION_H_
+#define STORM_STORAGE_BISIMULATION_PARTITION_H_
+
+#include <cstddef>
+#include <list>
+#include <memory>
+
+#include "src/storage/bisimulation/Block.h"
+
+#include "src/storage/BitVector.h"
+
+namespace storm {
+    namespace storage {
+        namespace bisimulation {
+            
+            template<typename DataType>
+            class Partition {
+            public:
+                /*!
+                 * Creates a partition with one block consisting of all the states.
+                 *
+                 * @param numberOfStates The number of states the partition holds.
+                 */
+                Partition(std::size_t numberOfStates);
+                
+                /*!
+                 * Creates a partition with three blocks: one with all phi states, one with all psi states and one with
+                 * all other states. The former two blocks are marked as being absorbing, because their outgoing
+                 * transitions shall not be taken into account for future refinement.
+                 *
+                 * @param numberOfStates The number of states the partition holds.
+                 * @param prob0States The states which have probability 0 of satisfying phi until psi.
+                 * @param prob1States The states which have probability 1 of satisfying phi until psi.
+                 * @param representativeProb1State If the set of prob1States is non-empty, this needs to be a state
+                 * that is representative for this block in the sense that the state representing this block in the quotient
+                 * model will receive exactly the atomic propositions of the representative state.
+                 */
+                Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, boost::optional<storm::storage::sparse::state_type> representativeProb1State);
+                
+                Partition() = default;
+                Partition(Partition const& other) = default;
+                Partition& operator=(Partition const& other) = default;
+                Partition(Partition&& other) = default;
+                Partition& operator=(Partition&& other) = default;
+                
+                // Retrieves the size of the partition, i.e. the number of blocks.
+                std::size_t size() const;
+                
+                // Prints the partition to the standard output.
+                void print() const;
+                
+                // Splits the block at the given position and inserts a new block after the current one holding the rest
+                // of the states.
+                std::pair<typename std::vector<std::unique_ptr<Block<DataType>>>::iterator, bool> splitBlock(Block<DataType>& block, storm::storage::sparse::state_type position);
+
+                // Sorts the given range of the partitition according to the given order.
+                void sortRange(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions = true);
+                
+                // Sorts the block according to the given order.
+                void sortBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions = true);
+
+                // Computes the start indices of equal ranges within the given range wrt. to the given less function.
+                std::vector<uint_fast64_t> computeRangesOfEqualValue(uint_fast64_t startIndex, uint_fast64_t endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less);
+                
+                // Splits the block by sorting the states according to the given function and then identifying the split
+                // points. The callback function is called for every newly created block.
+                bool splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback);
+
+                // Splits the block by sorting the states according to the given function and then identifying the split
+                // points.
+                bool splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less);
+                
+                // Splits all blocks by using the sorting-based splitting. The callback is called for all newly created
+                // blocks.
+                bool split(std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback);
+
+                // Splits all blocks by using the sorting-based splitting.
+                bool split(std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less);
+                
+                // Splits the block such that the resulting blocks contain only states in the given set or none at all.
+                // If the block is split, the given block will contain the states *not* in the given set and the newly
+                // created block will contain the states *in* the given set.
+                void splitStates(Block<DataType>& block, storm::storage::BitVector const& states);
+                
+                /*!
+                 * Splits all blocks of the partition such that afterwards all blocks contain only states within the given
+                 * set of states or no such state at all.
+                 */
+                void splitStates(storm::storage::BitVector const& states);
+                
+                // Sorts the block based on the state indices.
+                void sortBlock(Block<DataType> const& block);
+                
+                // Retrieves the blocks of the partition.
+                std::vector<std::unique_ptr<Block<DataType>>> const& getBlocks() const;
+                
+                // Retrieves the blocks of the partition.
+                std::vector<std::unique_ptr<Block<DataType>>>& getBlocks();
+                
+                // Checks the partition for internal consistency.
+                bool check() const;
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::iterator begin(Block<DataType> const& block);
+
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::const_iterator begin(Block<DataType> const& block) const;
+
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::iterator end(Block<DataType> const& block);
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::const_iterator end(Block<DataType> const& block) const;
+                
+                // Returns an iterator to the beginning of the states in the partition.
+                std::vector<storm::storage::sparse::state_type>::iterator begin();
+
+                // Returns an iterator to the beginning of the states in the partition.
+                std::vector<storm::storage::sparse::state_type>::const_iterator begin() const;
+
+                // Returns an iterator to the end of the states in the partition.
+                std::vector<storm::storage::sparse::state_type>::iterator end();
+                
+                // Returns an iterator to the end of the states in the partition.
+                std::vector<storm::storage::sparse::state_type>::const_iterator end() const;
+                
+                // Swaps the positions of the two given states.
+                void swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2);
+                
+                // Retrieves the block of the given state.
+                Block<DataType>& getBlock(storm::storage::sparse::state_type state);
+                
+                // Retrieves the block of the given state.
+                Block<DataType> const& getBlock(storm::storage::sparse::state_type state) const;
+                
+                // Retrieves the position of the given state.
+                storm::storage::sparse::state_type const& getPosition(storm::storage::sparse::state_type state) const;
+                
+                // Sets the position of the state to the given position.
+                storm::storage::sparse::state_type const& getState(storm::storage::sparse::state_type position) const;
+                
+                // Updates the block mapping for the given range of states to the specified block.
+                void mapStatesToBlock(Block<DataType>& block, std::vector<storm::storage::sparse::state_type>::iterator first, std::vector<storm::storage::sparse::state_type>::iterator last);
+                
+                // Update the state to position for the states in the given block.
+                void mapStatesToPositions(Block<DataType> const& block);
+
+                // Update the state to position for the states in the given range.
+                void mapStatesToPositions(std::vector<storm::storage::sparse::state_type>::const_iterator first, std::vector<storm::storage::sparse::state_type>::const_iterator last);
+                
+                // Swaps the positions of the two states given by their positions.
+                void swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2);
+                
+            private:
+                // FIXME: necessary?
+                // Inserts a block before the given block. The new block will cover all states between the beginning
+                // of the given block and the end of the previous block.
+                Block<DataType>& insertBlock(Block<DataType>& block);
+                
+                // FIXME: necessary?
+                // Sets the position of the given state.
+                void setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position);
+                
+                // The of blocks in the partition.
+                std::vector<std::unique_ptr<Block<DataType>>> blocks;
+                
+                // A mapping of states to their blocks.
+                std::vector<Block<DataType>*> stateToBlockMapping;
+                
+                // A vector containing all the states. It is ordered in a way such that the blocks only need to define
+                // their start/end indices.
+                std::vector<storm::storage::sparse::state_type> states;
+                
+                // This vector keeps track of the position of each state in the state vector.
+                std::vector<storm::storage::sparse::state_type> positions;
+            };
+        }
+    }
+}
+
+#endif /* STORM_STORAGE_BISIMULATION_PARTITION_H_ */
diff --git a/src/storage/expressions/Expression.cpp b/src/storage/expressions/Expression.cpp
index 725a41392..96cf57a31 100644
--- a/src/storage/expressions/Expression.cpp
+++ b/src/storage/expressions/Expression.cpp
@@ -162,6 +162,12 @@ namespace storm {
             return this->getBaseExpression().accept(visitor);
         }
         
+        std::string Expression::toString() {
+            std::stringstream stream;
+            stream << *this;
+            return stream.str();
+        }
+        
         std::ostream& operator<<(std::ostream& stream, Expression const& expression) {
             stream << expression.getBaseExpression();
             return stream;
diff --git a/src/storage/expressions/Expression.h b/src/storage/expressions/Expression.h
index c6ecdcbec..e8fe0e60b 100644
--- a/src/storage/expressions/Expression.h
+++ b/src/storage/expressions/Expression.h
@@ -292,6 +292,13 @@ namespace storm {
              */
             boost::any accept(ExpressionVisitor& visitor) const;
             
+            /*!
+             * Converts the expression into a string.
+             *
+             * @return The string representation of the expression.
+             */
+            std::string toString();
+            
             friend std::ostream& operator<<(std::ostream& stream, Expression const& expression);
 
         private:
diff --git a/src/storage/expressions/ToRationalFunctionVisitor.cpp b/src/storage/expressions/ToRationalFunctionVisitor.cpp
index 5740cc66d..97501aa1b 100644
--- a/src/storage/expressions/ToRationalFunctionVisitor.cpp
+++ b/src/storage/expressions/ToRationalFunctionVisitor.cpp
@@ -65,7 +65,7 @@ namespace storm {
             if (variablePair != variableToVariableMap.end()) {
                 return convertVariableToPolynomial(variablePair->second);
             } else {
-                carl::Variable carlVariable = carl::VariablePool::getInstance().getFreshVariable(expression.getVariableName());
+                carl::Variable carlVariable = carl::freshRealVariable(expression.getVariableName());
                 variableToVariableMap.emplace(expression.getVariable(), carlVariable);
                 return convertVariableToPolynomial(carlVariable);
             }
diff --git a/src/utility/ConstantsComparator.cpp b/src/utility/ConstantsComparator.cpp
index ed3d82fe4..4ad0a2612 100644
--- a/src/utility/ConstantsComparator.cpp
+++ b/src/utility/ConstantsComparator.cpp
@@ -35,6 +35,11 @@ namespace storm {
             return false;
         }
         
+        template<typename ValueType>
+        bool ConstantsComparator<ValueType>::isLess(ValueType const& value1, ValueType const& value2) const {
+            return value1 < value2;
+        }
+        
         ConstantsComparator<float>::ConstantsComparator() : precision(static_cast<float>(storm::settings::generalSettings().getPrecision())) {
             // Intentionally left empty.
         }
@@ -63,6 +68,10 @@ namespace storm {
             return value == storm::utility::infinity<float>();
         }
         
+        bool ConstantsComparator<float>::isLess(float const& value1, float const& value2) const {
+            return std::abs(value1 - value2) < precision;
+        }
+        
         ConstantsComparator<double>::ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) {
             // Intentionally left empty.
         }
@@ -90,7 +99,11 @@ namespace storm {
         bool ConstantsComparator<double>::isConstant(double const& value) const {
             return true;
         }
-                
+        
+        bool ConstantsComparator<double>::isLess(double const& value1, double const& value2) const {
+            return value1 < value2 - precision;
+        }
+        
         // Explicit instantiations.
         template class ConstantsComparator<double>;
         template class ConstantsComparator<float>;
diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h
index 7cd03eab3..d30a811d2 100644
--- a/src/utility/ConstantsComparator.h
+++ b/src/utility/ConstantsComparator.h
@@ -21,6 +21,8 @@ namespace storm {
             bool isConstant(ValueType const& value) const;
             
             bool isInfinity(ValueType const& value) const;
+            
+            bool isLess(ValueType const& value1, ValueType const& value2) const;
         };
         
         // For floats we specialize this class and consider the comparison modulo some predefined precision.
@@ -41,6 +43,8 @@ namespace storm {
             
             bool isInfinity(float const& value) const;
 
+            bool isLess(float const& value1, float const& value2) const;
+            
         private:
             // The precision used for comparisons.
             float precision;
@@ -64,6 +68,8 @@ namespace storm {
             
             bool isConstant(double const& value) const;
             
+            bool isLess(double const& value1, double const& value2) const;
+            
         private:
             // The precision used for comparisons.
             double precision;
diff --git a/src/utility/macros.h b/src/utility/macros.h
index b84cd53b5..037d1882b 100644
--- a/src/utility/macros.h
+++ b/src/utility/macros.h
@@ -12,19 +12,19 @@ extern log4cplus::Logger logger;
  * Define the macros STORM_LOG_DEBUG and STORM_LOG_TRACE.
  */
 #define STORM_LOG_DEBUG(message)                \
-{                                               \
+do {                                            \
     LOG4CPLUS_DEBUG(logger, message);           \
 } while (false)                                 \
 
 #define STORM_LOG_TRACE(message)                \
-{                                               \
+do {                                            \
     LOG4CPLUS_TRACE(logger, message);           \
 } while (false)                                 \
 
 // Define STORM_LOG_ASSERT which is only checked when NDEBUG is not set.
 #ifndef NDEBUG
 #define STORM_LOG_ASSERT(cond, message)         \
-{                                               \
+do {                                            \
 if (!(cond)) {                                  \
 LOG4CPLUS_ERROR(logger, message);               \
 assert(cond);                                   \
@@ -36,7 +36,7 @@ assert(cond);                                   \
 #endif
 // Define STORM_LOG_THROW to always throw the exception with the given message if the condition fails to hold.
 #define STORM_LOG_THROW(cond, exception, message)     \
-{                                               \
+do {                                            \
     if (!(cond)) {                              \
         LOG4CPLUS_ERROR(logger, message);       \
         throw exception() << message;           \
@@ -46,36 +46,36 @@ assert(cond);                                   \
 
 // Define STORM_LOG_WARN, STORM_LOG_ERROR and STORM_LOG_INFO to log the given message with the corresponding log levels.
 #define STORM_LOG_WARN(message)                 \
-{                                               \
+do {                                            \
     LOG4CPLUS_WARN(logger, message);            \
 } while (false)                                 \
 
 #define STORM_LOG_WARN_COND(cond, message)      \
-{                                               \
+do {                                            \
     if (!(cond)) {                              \
         LOG4CPLUS_WARN(logger, message);        \
     }                                           \
 } while (false)                                 \
 
 #define STORM_LOG_INFO(message)                 \
-{                                               \
+do {                                            \
     LOG4CPLUS_INFO(logger, message);            \
 } while (false)                                 \
 
 #define STORM_LOG_INFO_COND(cond, message)      \
-{                                               \
+do {                                            \
     if (!(cond)) {                              \
         LOG4CPLUS_INFO(logger, message);        \
     }                                           \
 } while (false)                                 \
 
 #define STORM_LOG_ERROR(message)                \
-{                                               \
+do {                                            \
     LOG4CPLUS_ERROR(logger, message);           \
 } while (false)                                 \
 
 #define STORM_LOG_ERROR_COND(cond, message)     \
-{                                               \
+do {                                            \
     if (!(cond)) {                              \
         LOG4CPLUS_ERROR(logger, message);       \
     }                                           \
diff --git a/src/utility/storm.h b/src/utility/storm.h
index b53647522..4a5d30f0b 100644
--- a/src/utility/storm.h
+++ b/src/utility/storm.h
@@ -30,6 +30,7 @@
 #include "src/models/ModelBase.h"
 #include "src/models/sparse/Model.h"
 #include "src/models/sparse/StandardRewardModel.h"
+#include "src/models/sparse/MarkovAutomaton.h"
 #include "src/models/symbolic/Model.h"
 #include "src/models/symbolic/StandardRewardModel.h"
 
@@ -43,7 +44,8 @@
 #include "src/builder/DdPrismModelBuilder.h"
 
 // Headers for model processing.
-#include "src/storage/DeterministicModelBisimulationDecomposition.h"
+#include "src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h"
+#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
 
 // Headers for model checking.
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
@@ -80,8 +82,6 @@ namespace storm {
     std::vector<std::shared_ptr<storm::logic::Formula>> parseFormulasForExplicit(std::string const& inputString);
     std::vector<std::shared_ptr<storm::logic::Formula>> parseFormulasForProgram(std::string const& inputString, storm::prism::Program const& program);
             
- 
- 
     template<typename ValueType>
     std::shared_ptr<storm::models::ModelBase> buildSymbolicModel(storm::prism::Program const& program, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
         std::shared_ptr<storm::models::ModelBase> result(nullptr);
@@ -111,37 +111,74 @@ namespace storm {
             result = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
         }
 
-        // Then, build the model from the symbolic description.
         return result;
     }
-            
-    template<typename ValueType>
-    std::shared_ptr<storm::models::ModelBase> preprocessModel(std::shared_ptr<storm::models::ModelBase> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+    
+    template<typename ModelType>
+    std::shared_ptr<ModelType> performDeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+        std::cout << "Performing bisimulation minimization... ";
+        typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options options;
+        if (!formulas.empty()) {
+            options = typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options(*model, formulas);
+        }
+        if (storm::settings::bisimulationSettings().isWeakBisimulationSet()) {
+            options.type = storm::storage::BisimulationType::Weak;
+            options.bounded = false;
+        }
+        
+        storm::storage::DeterministicModelBisimulationDecomposition<ModelType> bisimulationDecomposition(*model, options);
+        bisimulationDecomposition.computeBisimulationDecomposition();
+        model = bisimulationDecomposition.getQuotient();
+        std::cout << "done." << std::endl << std::endl;
+        return model;
+    }
+    
+    template<typename ModelType>
+    std::shared_ptr<ModelType> performNondeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+        std::cout << "Performing bisimulation minimization... ";
+        typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options options;
+        if (!formulas.empty()) {
+            options = typename storm::storage::NondeterministicModelBisimulationDecomposition<ModelType>::Options(*model, formulas);
+        }
+        if (storm::settings::bisimulationSettings().isWeakBisimulationSet()) {
+            options.type = storm::storage::BisimulationType::Weak;
+            options.bounded = false;
+        }
+        
+        storm::storage::NondeterministicModelBisimulationDecomposition<ModelType> bisimulationDecomposition(*model, options);
+        bisimulationDecomposition.computeBisimulationDecomposition();
+        model = bisimulationDecomposition.getQuotient();
+        std::cout << "done." << std::endl << std::endl;
+        return model;
+    }
+    
+    template<typename ModelType, typename ValueType = typename ModelType::ValueType, typename std::enable_if<std::is_base_of<storm::models::sparse::Model<ValueType>, ModelType>::value, bool>::type = 0>
+    std::shared_ptr<storm::models::ModelBase> preprocessModel(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
         if (storm::settings::generalSettings().isBisimulationSet()) {
+            
             STORM_LOG_THROW(model->isSparseModel(), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for sparse models.");
-            std::shared_ptr<storm::models::sparse::Model<ValueType>> sparseModel = model->template as<storm::models::sparse::Model<ValueType>>();
-            STORM_LOG_THROW(model->getType() == storm::models::ModelType::Dtmc || model->getType() == storm::models::ModelType::Ctmc, storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs.");
-            std::shared_ptr<storm::models::sparse::Dtmc<ValueType>> dtmc = sparseModel->template as<storm::models::sparse::Dtmc<ValueType>>();
+            STORM_LOG_THROW(model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc) || model->isOfType(storm::models::ModelType::Mdp), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs, CTMCs and MDPs.");
 
-            dtmc->reduceToStateBasedRewards();
+            model->reduceToStateBasedRewards();
 
-            std::cout << "Performing bisimulation minimization... ";
-            typename storm::storage::DeterministicModelBisimulationDecomposition<ValueType>::Options options;
-            if (!formulas.empty()) {
-                options = typename storm::storage::DeterministicModelBisimulationDecomposition<ValueType>::Options(*sparseModel, formulas);
-            }
-            if (storm::settings::bisimulationSettings().isWeakBisimulationSet()) {
-                options.weak = true;
-                options.bounded = false;
+            if (model->isOfType(storm::models::ModelType::Dtmc)) {
+                return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Dtmc<ValueType>>(model->template as<storm::models::sparse::Dtmc<ValueType>>(), formulas);
+            } else if (model->isOfType(storm::models::ModelType::Ctmc)) {
+                return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas);
+            } else {
+                return performNondeterministicSparseBisimulationMinimization<storm::models::sparse::Mdp<ValueType>>(model->template as<storm::models::sparse::Mdp<ValueType>>(), formulas);
             }
-
-            storm::storage::DeterministicModelBisimulationDecomposition<ValueType> bisimulationDecomposition(*dtmc, options);
-            model = bisimulationDecomposition.getQuotient();
-            std::cout << "done." << std::endl << std::endl;
+            
         }
         return model;
     }
-            
+    
+    template<typename ModelType, storm::dd::DdType DdType = ModelType::DdType, typename std::enable_if< std::is_base_of< storm::models::symbolic::Model<DdType>, ModelType >::value, bool>::type = 0>
+    std::shared_ptr<storm::models::ModelBase> preprocessModel(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+        // No preprocessing available yet.
+        return model;
+    }
+    
     template<typename ValueType>
     void generateCounterexample(storm::prism::Program const& program, std::shared_ptr<storm::models::sparse::Model<ValueType>> model, std::shared_ptr<storm::logic::Formula> const& formula) {
         if (storm::settings::counterexampleGeneratorSettings().isMinimalCommandSetGenerationSet()) {
@@ -231,7 +268,6 @@ namespace storm {
         return result;
 
     }
-
             
 #ifdef STORM_HAVE_CARL
     inline void exportParametricResultToFile(storm::RationalFunction const& result, storm::models::sparse::Dtmc<storm::RationalFunction>::ConstraintCollector const& constraintCollector, std::string const& path) {
@@ -264,7 +300,6 @@ namespace storm {
         return result;
     }
             
-
 #endif
 
     template<storm::dd::DdType DdType>
diff --git a/test/functional/storage/DeterministicModelBisimulationDecompositionTest.cpp b/test/functional/storage/DeterministicModelBisimulationDecompositionTest.cpp
index 5081e9e02..a706475bd 100644
--- a/test/functional/storage/DeterministicModelBisimulationDecompositionTest.cpp
+++ b/test/functional/storage/DeterministicModelBisimulationDecompositionTest.cpp
@@ -1,7 +1,8 @@
 #include "gtest/gtest.h"
 #include "storm-config.h"
 #include "src/parser/AutoParser.h"
-#include "src/storage/DeterministicModelBisimulationDecomposition.h"
+#include "src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h"
+#include "src/models/sparse/Dtmc.h"
 #include "src/models/sparse/StandardRewardModel.h"
 
 TEST(DeterministicModelBisimulationDecomposition, Die) {
@@ -10,7 +11,8 @@ TEST(DeterministicModelBisimulationDecomposition, Die) {
     ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc);
     std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = abstractModel->as<storm::models::sparse::Dtmc<double>>();
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim(*dtmc);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim(*dtmc);
+    ASSERT_NO_THROW(bisim.computeBisimulationDecomposition());
     std::shared_ptr<storm::models::sparse::Model<double>> result;
     ASSERT_NO_THROW(result = bisim.getQuotient());
 
@@ -19,13 +21,14 @@ TEST(DeterministicModelBisimulationDecomposition, Die) {
     EXPECT_EQ(20ul, result->getNumberOfTransitions());
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options;
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options;
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options;
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options;
 #endif
     options.respectedAtomicPropositions = std::set<std::string>({"one"});
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim2(*dtmc, options);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim2(*dtmc, options);
+    ASSERT_NO_THROW(bisim2.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim2.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -33,9 +36,10 @@ TEST(DeterministicModelBisimulationDecomposition, Die) {
     EXPECT_EQ(8ul, result->getNumberOfTransitions());
 
     options.bounded = false;
-    options.weak = true;
+    options.type = storm::storage::BisimulationType::Weak;
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim3(*dtmc, options);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim3(*dtmc, options);
+    ASSERT_NO_THROW(bisim3.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim3.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -46,12 +50,13 @@ TEST(DeterministicModelBisimulationDecomposition, Die) {
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options2(*dtmc, *eventuallyFormula);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options2(*dtmc, *eventuallyFormula);
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options2(*dtmc, *eventuallyFormula);
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options2(*dtmc, *eventuallyFormula);
 #endif
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim4(*dtmc, options2);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim4(*dtmc, options2);
+    ASSERT_NO_THROW(bisim4.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim4.getQuotient());
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
     EXPECT_EQ(5ul, result->getNumberOfStates());
@@ -64,8 +69,9 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc);
     std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = abstractModel->as<storm::models::sparse::Dtmc<double>>();
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim(*dtmc);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim(*dtmc);
     std::shared_ptr<storm::models::sparse::Model<double>> result;
+    ASSERT_NO_THROW(bisim.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -73,13 +79,14 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     EXPECT_EQ(546ul, result->getNumberOfTransitions());
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options;
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options;
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options;
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options;
 #endif
     options.respectedAtomicPropositions = std::set<std::string>({"observe0Greater1"});
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim2(*dtmc, options);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim2(*dtmc, options);
+    ASSERT_NO_THROW(bisim2.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim2.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -87,9 +94,10 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     EXPECT_EQ(105ul, result->getNumberOfTransitions());
 
     options.bounded = false;
-    options.weak = true;
+    options.type = storm::storage::BisimulationType::Weak;
 
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim3(*dtmc, options);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim3(*dtmc, options);
+    ASSERT_NO_THROW(bisim3.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim3.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -100,11 +108,12 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options2(*dtmc, *eventuallyFormula);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options2(*dtmc, *eventuallyFormula);
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options2(*dtmc, *eventuallyFormula);
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options2(*dtmc, *eventuallyFormula);
 #endif
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim4(*dtmc, options2);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim4(*dtmc, options2);
+    ASSERT_NO_THROW(bisim4.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim4.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -114,11 +123,12 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     auto probabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(eventuallyFormula);
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options3(*dtmc, *probabilityOperatorFormula);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options3(*dtmc, *probabilityOperatorFormula);
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options3(*dtmc, *probabilityOperatorFormula);
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options3(*dtmc, *probabilityOperatorFormula);
 #endif
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim5(*dtmc, options3);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim5(*dtmc, options3);
+    ASSERT_NO_THROW(bisim5.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim5.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
@@ -128,11 +138,12 @@ TEST(DeterministicModelBisimulationDecomposition, Crowds) {
     auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(std::make_shared<storm::logic::BooleanLiteralFormula>(true), labelFormula, 50);
 
 #ifdef WINDOWS
-    storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options4(*dtmc, *boundedUntilFormula);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options4(*dtmc, *boundedUntilFormula);
 #else
-    typename storm::storage::DeterministicModelBisimulationDecomposition<double>::Options options4(*dtmc, *boundedUntilFormula);
+    typename storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>>::Options options4(*dtmc, *boundedUntilFormula);
 #endif
-    storm::storage::DeterministicModelBisimulationDecomposition<double> bisim6(*dtmc, options4);
+    storm::storage::DeterministicModelBisimulationDecomposition<storm::models::sparse::Dtmc<double>> bisim6(*dtmc, options4);
+    ASSERT_NO_THROW(bisim6.computeBisimulationDecomposition());
     ASSERT_NO_THROW(result = bisim6.getQuotient());
 
     EXPECT_EQ(storm::models::ModelType::Dtmc, result->getType());
diff --git a/test/functional/storage/NondeterministicModelBisimulationDecompositionTest.cpp b/test/functional/storage/NondeterministicModelBisimulationDecompositionTest.cpp
new file mode 100644
index 000000000..e2207bebb
--- /dev/null
+++ b/test/functional/storage/NondeterministicModelBisimulationDecompositionTest.cpp
@@ -0,0 +1,58 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/parser/PrismParser.h"
+#include "src/parser/FormulaParser.h"
+
+#include "src/builder/ExplicitPrismModelBuilder.h"
+
+#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
+#include "src/models/sparse/Mdp.h"
+#include "src/models/sparse/StandardRewardModel.h"
+
+TEST(NondeterministicModelBisimulationDecomposition, TwoDice) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/two_dice.nm");
+
+    // Build the die model without its reward model.
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    std::shared_ptr<storm::models::sparse::Mdp<double>> mdp = model->as<storm::models::sparse::Mdp<double>>();
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim(*mdp);
+    ASSERT_NO_THROW(bisim.computeBisimulationDecomposition());
+    std::shared_ptr<storm::models::sparse::Model<double>> result;
+    ASSERT_NO_THROW(result = bisim.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(77ul, result->getNumberOfStates());
+    EXPECT_EQ(183ul, result->getNumberOfTransitions());
+    EXPECT_EQ(97ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options;
+    options.respectedAtomicPropositions = std::set<std::string>({"two"});
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim2(*mdp, options);
+    ASSERT_NO_THROW(bisim2.computeBisimulationDecomposition());
+    ASSERT_NO_THROW(result = bisim2.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(11ul, result->getNumberOfStates());
+    EXPECT_EQ(26ul, result->getNumberOfTransitions());
+    EXPECT_EQ(14ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    std::shared_ptr<storm::logic::Formula> formula = formulaParser.parseSingleFormulaFromString("Pmin=? [F \"two\"]");
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options2(*mdp, *formula);
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim3(*mdp, options2);
+    ASSERT_NO_THROW(bisim3.computeBisimulationDecomposition());
+    ASSERT_NO_THROW(result = bisim3.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(11ul, result->getNumberOfStates());
+    EXPECT_EQ(26ul, result->getNumberOfTransitions());
+    EXPECT_EQ(14ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+}
\ No newline at end of file