diff --git a/src/models/AtomicPropositionsLabeling.h b/src/models/AtomicPropositionsLabeling.h
index 1969a525a..8959756a4 100644
--- a/src/models/AtomicPropositionsLabeling.h
+++ b/src/models/AtomicPropositionsLabeling.h
@@ -250,6 +250,19 @@ public:
 		return this->singleLabelings[apIndexPair->second].get(state);
 	}
 
+    /*!
+     * Retrieves the set of atomic propositions of the model.
+     *
+     * @return The set of atomic propositions of the model.
+     */
+    std::set<std::string> getAtomicPropositions() const {
+        std::set<std::string> result;
+        for (auto const& labeling : this->nameToLabelingMap) {
+            result.insert(labeling.first);
+        }
+        return result;
+    }
+    
 	/*!
 	 * Returns the number of atomic propositions managed by this object (set in the initialization).
      *
diff --git a/src/models/Ctmc.h b/src/models/Ctmc.h
index 229ccfe9e..b3c1f1ef1 100644
--- a/src/models/Ctmc.h
+++ b/src/models/Ctmc.h
@@ -36,8 +36,9 @@ public:
 	 * propositions to each state.
 	 */
 	Ctmc(storm::storage::SparseMatrix<T> const& rateMatrix, storm::models::AtomicPropositionsLabeling const& stateLabeling,
-				boost::optional<std::vector<T>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>> const& optionalTransitionRewardMatrix,
-                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling)
+				boost::optional<std::vector<T>> const& optionalStateRewardVector = boost::optional<std::vector<T>>(),
+                boost::optional<storm::storage::SparseMatrix<T>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<T>>(),
+                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>())
 			: AbstractDeterministicModel<T>(rateMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
 		// Intentionally left empty.
 	}
@@ -51,8 +52,9 @@ public:
 	 * propositions to each state.
 	 */
 	Ctmc(storm::storage::SparseMatrix<T>&& rateMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling,
-				boost::optional<std::vector<T>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,
-                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+				boost::optional<std::vector<T>>&& optionalStateRewardVector = boost::optional<std::vector<T>>(),
+                boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<T>>(),
+                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>())
 			// The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class
 			: AbstractDeterministicModel<T>(std::move(rateMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix),
                                             std::move(optionalChoiceLabeling)) {
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index dcb804fe8..456d3a179 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -73,8 +73,9 @@ public:
 	 * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
 	 */
 	Dtmc(storm::storage::SparseMatrix<T>&& probabilityMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling,
-				boost::optional<std::vector<T>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,
-                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling)
+				boost::optional<std::vector<T>>&& optionalStateRewardVector = boost::optional<std::vector<T>>(),
+                boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<T>>(),
+                boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>())
 				// The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class
 			: AbstractDeterministicModel<T>(std::move(probabilityMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix),
                                             std::move(optionalChoiceLabeling)) {
diff --git a/src/settings/modules/GeneralSettings.cpp b/src/settings/modules/GeneralSettings.cpp
index 0e1b8775a..d04ec0d25 100644
--- a/src/settings/modules/GeneralSettings.cpp
+++ b/src/settings/modules/GeneralSettings.cpp
@@ -43,6 +43,8 @@ namespace storm {
             const std::string GeneralSettings::constantsOptionShortName = "const";
             const std::string GeneralSettings::statisticsOptionName = "statistics";
             const std::string GeneralSettings::statisticsOptionShortName = "stats";
+            const std::string GeneralSettings::bisimulationOptionName = "bisimulation";
+            const std::string GeneralSettings::bisimulationOptionShortName = "bisim";
 
             GeneralSettings::GeneralSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) {
                 this->addOption(storm::settings::OptionBuilder(moduleName, helpOptionName, false, "Shows all available options, arguments and descriptions.").setShortName(helpOptionShortName)
@@ -74,7 +76,7 @@ namespace storm {
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the LTL formulas.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, counterexampleOptionName, false, "Generates a counterexample for the given PRCTL formulas if not satisfied by the model")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file to which the counterexample is to be written.").setDefaultValueString("-").setIsOptional(true).build()).setShortName(counterexampleOptionShortName).build());
-                
+                this->addOption(storm::settings::OptionBuilder(moduleName, bisimulationOptionName, false, "Sets whether to perform bisimulation minimization.").setShortName(bisimulationOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, transitionRewardsOptionName, "", "If given, the transition rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the transition rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, stateRewardsOptionName, false, "If given, the state rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
@@ -93,7 +95,6 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, constantsOptionName, false, "Specifies the constant replacements to use in symbolic models. Note that Note that this requires the model to be given as an symbolic model (i.e., via --" + symbolicOptionName + ").").setShortName(constantsOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of constants and their value, e.g. a=1,b=2,c=3.").setDefaultValueString("").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, statisticsOptionName, false, "Sets whether to display statistics if available.").setShortName(statisticsOptionShortName).build());
-
             }
             
             bool GeneralSettings::isHelpSet() const {
@@ -284,6 +285,10 @@ namespace storm {
                 return true;
             }
             
+            bool GeneralSettings::isBisimulationSet() const {
+                return this->getOption(bisimulationOptionName).getHasOptionBeenSet();
+            }
+            
         } // namespace modules
     } // namespace settings
 } // namespace storm
\ No newline at end of file
diff --git a/src/settings/modules/GeneralSettings.h b/src/settings/modules/GeneralSettings.h
index 912053083..c65ab57a2 100644
--- a/src/settings/modules/GeneralSettings.h
+++ b/src/settings/modules/GeneralSettings.h
@@ -322,6 +322,13 @@ namespace storm {
                  */
                 bool isShowStatisticsSet() const;
                 
+                /*!
+                 * Retrieves whether the option to perform bisimulation minimization is set.
+                 *
+                 * @return True iff the option was set.
+                 */
+                bool isBisimulationSet() const;
+                
                 bool check() const override;
 
                 // The name of the module.
@@ -363,6 +370,8 @@ namespace storm {
                 static const std::string constantsOptionShortName;
                 static const std::string statisticsOptionName;
                 static const std::string statisticsOptionShortName;
+                static const std::string bisimulationOptionName;
+                static const std::string bisimulationOptionShortName;
             };
             
         } // namespace modules
diff --git a/src/storage/Decomposition.cpp b/src/storage/Decomposition.cpp
index 06656b483..39c49bf4f 100644
--- a/src/storage/Decomposition.cpp
+++ b/src/storage/Decomposition.cpp
@@ -89,7 +89,10 @@ namespace storm {
             out << "]";
             return out;
         }
-        
+
+        template class Decomposition<StateBlock>;
+        template std::ostream& operator<<(std::ostream& out, Decomposition<StateBlock> const& decomposition);
+
         template class Decomposition<StronglyConnectedComponent>;
         template std::ostream& operator<<(std::ostream& out, Decomposition<StronglyConnectedComponent> const& decomposition);
         
diff --git a/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp
new file mode 100644
index 000000000..93433d19f
--- /dev/null
+++ b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp
@@ -0,0 +1,904 @@
+#include "src/storage/DeterministicModelStrongBisimulationDecomposition.h"
+
+#include <algorithm>
+#include <unordered_map>
+#include <chrono>
+#include <iomanip>
+
+#include "src/utility/graph.h"
+#include "src/exceptions/IllegalFunctionCallException.h"
+
+namespace storm {
+    namespace storage {
+        
+        template<typename ValueType>
+        std::size_t DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::blockId = 0;
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::shared_ptr<std::string> const& label) : next(next), prev(prev), begin(begin), end(end), markedAsSplitter(false), markedAsPredecessorBlock(false), markedPosition(begin), absorbing(false), id(blockId++), label(label) {
+            if (next != nullptr) {
+                next->prev = this;
+            }
+            if (prev != nullptr) {
+                prev->next = this;
+            }
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<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.states[index] << " ";
+            }
+            std::cout << std::endl << "original: " << std::endl;
+            for (storm::storage::sparse::state_type index = this->getOriginalBegin(); index < this->end; ++index) {
+                std::cout << partition.states[index] << " ";
+            }
+            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.values[index] << " ";
+            }
+            std::cout << std::endl;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::setBegin(storm::storage::sparse::state_type begin) {
+            this->begin = begin;
+            this->markedPosition = begin;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::setEnd(storm::storage::sparse::state_type end) {
+            this->end = end;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::incrementBegin() {
+            ++this->begin;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::decrementEnd() {
+            ++this->begin;
+        }
+        
+        template<typename ValueType>
+        storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getBegin() const {
+            return this->begin;
+        }
+        
+        template<typename ValueType>
+        storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getOriginalBegin() const {
+            if (this->hasPreviousBlock()) {
+                return this->getPreviousBlock().getEnd();
+            } else {
+                return 0;
+            }
+        }
+        
+        template<typename ValueType>
+        storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getEnd() const {
+            return this->end;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::setIterator(const_iterator it) {
+            this->selfIt = it;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::const_iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getIterator() const {
+            return this->selfIt;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::const_iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getNextIterator() const {
+            return this->getNextBlock().getIterator();
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::const_iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getPreviousIterator() const {
+            return this->getPreviousBlock().getIterator();
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getNextBlock() {
+            return *this->next;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getNextBlock() const {
+            return *this->next;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::hasNextBlock() const {
+            return this->next != nullptr;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getPreviousBlock() {
+            return *this->prev;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block* DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getPreviousBlockPointer() {
+            return this->prev;
+        }
+
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block* DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getNextBlockPointer() {
+            return this->next;
+        }
+
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getPreviousBlock() const {
+            return *this->prev;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::hasPreviousBlock() const {
+            return this->prev != nullptr;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::check() const {
+            if (this->begin >= this->end) {
+                assert(false);
+            }
+            if (this->prev != nullptr) {
+                assert (this->prev->next == this);
+            }
+            if (this->next != nullptr) {
+                assert (this->next->prev == this);
+            }
+            return true;
+        }
+        
+        template<typename ValueType>
+        std::size_t DeterministicModelStrongBisimulationDecomposition<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 DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::isMarkedAsSplitter() const {
+            return this->markedAsSplitter;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::markAsSplitter() {
+            this->markedAsSplitter = true;
+        }
+
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::unmarkAsSplitter() {
+            this->markedAsSplitter = false;
+        }
+
+        template<typename ValueType>
+        std::size_t DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getId() const {
+            return this->id;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::setAbsorbing(bool absorbing) {
+            this->absorbing = absorbing;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::isAbsorbing() const {
+            return this->absorbing;
+        }
+        
+        template<typename ValueType>
+        storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getMarkedPosition() const {
+            return this->markedPosition;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::setMarkedPosition(storm::storage::sparse::state_type position) {
+            this->markedPosition = position;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::resetMarkedPosition() {
+            this->markedPosition = this->begin;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::incrementMarkedPosition() {
+            ++this->markedPosition;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::markAsPredecessorBlock() {
+            this->markedAsPredecessorBlock = true;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::unmarkAsPredecessorBlock() {
+            this->markedAsPredecessorBlock = false;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::isMarkedAsPredecessor() const {
+            return markedAsPredecessorBlock;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::hasLabel() const {
+            return this->label != nullptr;
+        }
+        
+        template<typename ValueType>
+        std::string const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getLabel() const {
+            STORM_LOG_THROW(this->label != nullptr, storm::exceptions::IllegalFunctionCallException, "Unable to retrieve label of block that has none.");
+            return *this->label;
+        }
+        
+        template<typename ValueType>
+        std::shared_ptr<std::string> const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Block::getLabelPtr() const {
+            return this->label;
+        }
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::Partition(std::size_t numberOfStates) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates), values(numberOfStates) {
+            // 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);
+            it->setIterator(it);
+            
+            // 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();
+            }
+        }
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, std::string const& otherLabel, std::string const& prob1Label) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates), values(numberOfStates) {
+            typename std::list<Block>::iterator firstIt = blocks.emplace(this->blocks.end(), 0, prob0States.getNumberOfSetBits(), nullptr, nullptr);
+            Block& firstBlock = *firstIt;
+            firstBlock.setIterator(firstIt);
+
+            storm::storage::sparse::state_type position = 0;
+            for (auto state : prob0States) {
+                states[position] = state;
+                positions[state] = position;
+                stateToBlockMapping[state] = &firstBlock;
+                ++position;
+            }
+            firstBlock.setAbsorbing(true);
+
+            typename std::list<Block>::iterator secondIt = blocks.emplace(this->blocks.end(), position, position + prob1States.getNumberOfSetBits(), &firstBlock, nullptr, std::shared_ptr<std::string>(new std::string(prob1Label)));
+            Block& secondBlock = *secondIt;
+            secondBlock.setIterator(secondIt);
+            
+            for (auto state : prob1States) {
+                states[position] = state;
+                positions[state] = position;
+                stateToBlockMapping[state] = &secondBlock;
+                ++position;
+            }
+            secondBlock.setAbsorbing(true);
+            
+            typename std::list<Block>::iterator thirdIt = blocks.emplace(this->blocks.end(), position, numberOfStates, &secondBlock, nullptr, otherLabel == "true" ? std::shared_ptr<std::string>(nullptr) : std::shared_ptr<std::string>(new std::string(otherLabel)));
+            Block& thirdBlock = *thirdIt;
+            thirdBlock.setIterator(thirdIt);
+            
+            storm::storage::BitVector otherStates = ~(prob0States | prob1States);
+            for (auto state : otherStates) {
+                states[position] = state;
+                positions[state] = position;
+                stateToBlockMapping[state] = &thirdBlock;
+                ++position;
+            }
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::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->values[this->positions[state1]], this->values[this->positions[state2]]);
+            std::swap(this->positions[state1], this->positions[state2]);
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::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]);
+            std::swap(this->values[position1], this->values[position2]);
+            
+            this->positions[state1] = position2;
+            this->positions[state2] = position1;
+        }
+        
+        template<typename ValueType>
+        storm::storage::sparse::state_type const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getPosition(storm::storage::sparse::state_type state) const {
+            return this->positions[state];
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<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& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getState(storm::storage::sparse::state_type position) const {
+            return this->states[position];
+        }
+        
+        template<typename ValueType>
+        ValueType const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getValue(storm::storage::sparse::state_type state) const {
+            return this->values[this->positions[state]];
+        }
+        
+        template<typename ValueType>
+        ValueType const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getValueAtPosition(storm::storage::sparse::state_type position) const {
+            return this->values[position];
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::setValue(storm::storage::sparse::state_type state, ValueType value) {
+            this->values[this->positions[state]] = value;
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType>& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getValues() {
+            return this->values;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::increaseValue(storm::storage::sparse::state_type state, ValueType value) {
+            this->values[this->positions[state]] += value;
+        }
+
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::updateBlockMapping(Block& 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 ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getFirstBlock() {
+            return *this->blocks.begin();
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBlock(storm::storage::sparse::state_type state) {
+            return *this->stateToBlockMapping[state];
+        }
+
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBlock(storm::storage::sparse::state_type state) const {
+            return *this->stateToBlockMapping[state];
+        }
+
+        template<typename ValueType>
+        std::vector<storm::storage::sparse::state_type>::iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBeginOfStates(Block const& block) {
+            return this->states.begin() + block.getBegin();
+        }
+        
+        template<typename ValueType>
+        std::vector<storm::storage::sparse::state_type>::iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getEndOfStates(Block const& block) {
+            return this->states.begin() + block.getEnd();
+        }
+
+        template<typename ValueType>
+        std::vector<storm::storage::sparse::state_type>::const_iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBeginOfStates(Block const& block) const {
+            return this->states.begin() + block.getBegin();
+        }
+        
+        template<typename ValueType>
+        std::vector<storm::storage::sparse::state_type>::const_iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getEndOfStates(Block const& block) const {
+            return this->states.begin() + block.getEnd();
+        }
+        
+        template<typename ValueType>
+        typename std::vector<ValueType>::iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBeginOfValues(Block const& block) {
+            return this->values.begin() + block.getBegin();
+        }
+        
+        template<typename ValueType>
+        typename std::vector<ValueType>::iterator DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getEndOfValues(Block const& block) {
+            return this->values.begin() + block.getEnd();
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::splitBlock(Block& block, storm::storage::sparse::state_type position) {
+            // FIXME: this could be improved by splitting off the smaller of the two resulting blocks.
+            
+            // 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, block.getLabelPtr());
+            selfIt->setIterator(selfIt);
+            Block& newBlock = *selfIt;
+            
+            // Make the current block end at the given position.
+            block.setBegin(position);
+            
+            // Update the mapping of the states in the newly created block.
+            for (auto it = this->getBeginOfStates(newBlock), ite = this->getEndOfStates(newBlock); it != ite; ++it) {
+                stateToBlockMapping[*it] = &newBlock;
+            }
+            
+            return newBlock;
+        }
+        
+        template<typename ValueType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block& DeterministicModelStrongBisimulationDecomposition<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);
+            Block& newBlock = *it;
+            newBlock.setIterator(it);
+            
+            // Update the mapping of the states in the newly created block.
+            for (auto it = this->getBeginOfStates(newBlock), ite = this->getEndOfStates(newBlock); it != ite; ++it) {
+                stateToBlockMapping[*it] = &newBlock;
+            }
+            
+            return *it;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<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->getBeginOfStates(block), this->getEndOfStates(block), [&statesWithLabel] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statesWithLabel.get(a) && !statesWithLabel.get(b); } );
+                
+                // Update the positions vector.
+                storm::storage::sparse::state_type position = block.getBegin();
+                for (auto stateIt = this->getBeginOfStates(block), stateIte = this->getEndOfStates(block); stateIt != stateIte; ++stateIt, ++position) {
+                    this->positions[*stateIt] = position;
+                }
+                
+                // Now we can find the first position in the block that does not have the label and create new blocks.
+                std::vector<storm::storage::sparse::state_type>::iterator it = std::find_if(this->getBeginOfStates(block), this->getEndOfStates(block), [&] (storm::storage::sparse::state_type const& a) { return !statesWithLabel.get(a); });
+                
+                // If not all the states agreed on the validity of the label, we need to split the block.
+                if (it != this->getBeginOfStates(block) && it != this->getEndOfStates(block)) {
+                    auto cutPoint = std::distance(this->states.begin(), it);
+                    this->splitBlock(block, cutPoint);
+                } else {
+                    // Otherwise, we simply proceed to the next block.
+                    ++blockIterator;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        std::list<typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block> const& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBlocks() const {
+            return this->blocks;
+        }
+        
+        template<typename ValueType>
+        std::vector<storm::storage::sparse::state_type>& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getStates() {
+            return this->states;
+        }
+        
+        template<typename ValueType>
+        std::list<typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Block>& DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::getBlocks() {
+            return this->blocks;
+        }
+        
+        template<typename ValueType>
+        bool DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::check() const {
+            for (uint_fast64_t state = 0; state < this->positions.size(); ++state) {
+                if (this->states[this->positions[state]] != state) {
+                    assert(false);
+                }
+            }
+            for (auto const& block : this->blocks) {
+                assert(block.check());
+                for (auto stateIt = this->getBeginOfStates(block), stateIte = this->getEndOfStates(block); stateIt != stateIte; ++stateIt) {
+                    if (this->stateToBlockMapping[*stateIt] != &block) {
+                        assert(false);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::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 << " ";
+            }
+            std::cout << std::endl;
+            std::cout << "size: " << this->blocks.size() << std::endl;
+            assert(this->check());
+        }
+        
+        template<typename ValueType>
+        std::size_t DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition::size() const {
+            return this->blocks.size();
+        }
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc<ValueType> const& model, bool buildQuotient) {
+            STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures.");
+            Partition initialPartition = getLabelBasedInitialPartition(model);
+            partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient);
+        }
+
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc<ValueType> const& model, bool buildQuotient) {
+            STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures.");
+            Partition initialPartition = getLabelBasedInitialPartition(model);
+            partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient);
+        }
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc<ValueType> const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient) {
+            STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures.");
+            storm::storage::SparseMatrix<ValueType> backwardTransitions = model.getBackwardTransitions();
+            Partition initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, phiLabel, psiLabel, bounded);
+            partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient);
+        }
+        
+        template<typename ValueType>
+        DeterministicModelStrongBisimulationDecomposition<ValueType>::DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc<ValueType> const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient) {
+            STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures.");
+            storm::storage::SparseMatrix<ValueType> backwardTransitions = model.getBackwardTransitions();
+            Partition initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, phiLabel, psiLabel, bounded);
+            partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient);
+        }
+        
+        template<typename ValueType>
+        std::shared_ptr<storm::models::AbstractDeterministicModel<ValueType>> DeterministicModelStrongBisimulationDecomposition<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 DeterministicModelStrongBisimulationDecomposition<ValueType>::buildQuotient(ModelType const& model, Partition const& partition) {
+            // 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::AtomicPropositionsLabeling newLabeling(this->size(), model.getStateLabeling().getNumberOfAtomicPropositions());
+            std::set<std::string> atomicPropositionsSet = model.getStateLabeling().getAtomicPropositions();
+            std::vector<std::string> atomicPropositions = std::vector<std::string>(atomicPropositionsSet.begin(), atomicPropositionsSet.end());
+            for (auto const& ap : atomicPropositions) {
+                newLabeling.addAtomicProposition(ap);
+            }
+            
+            // 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. It doesn't matter which state it is, because they all behave equally.
+                storm::storage::sparse::state_type representativeState = *block.begin();
+                
+                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::constantOne<ValueType>());
+                    
+                    if (oldBlock.hasLabel()) {
+                        newLabeling.addAtomicPropositionToState(oldBlock.getLabel(), 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();
+                        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) {
+                        builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second);
+                    }
+                    
+                    // If the block has a special label, we only add that to the block.
+                    if (oldBlock.hasLabel()) {
+                        newLabeling.addAtomicPropositionToState(oldBlock.getLabel(), blockIndex);
+                    } else {
+                        // Otherwise add all atomic propositions to the equivalence class that the representative state
+                        // satisfies.
+                        for (auto const& ap : atomicPropositions) {
+                            if (model.getStateLabeling().getStateHasAtomicProposition(ap, representativeState)) {
+                                newLabeling.addAtomicPropositionToState(ap, blockIndex);
+                            }
+                        }
+                    }
+                }
+            }
+            
+            // 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.addAtomicPropositionToState("init", initialBlock.getId());
+            }
+            
+            // FIXME:
+            // If reward structures are allowed, the quotient structures need to be built here.
+            
+            // Finally construct the quotient model.
+            this->quotient = std::shared_ptr<storm::models::AbstractDeterministicModel<ValueType>>(new ModelType(builder.build(), std::move(newLabeling)));
+        }
+        
+        template<typename ValueType>
+        template<typename ModelType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::partitionRefinement(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition, bool buildQuotient) {
+            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::deque<Block*> splitterQueue;
+            std::for_each(partition.getBlocks().begin(), partition.getBlocks().end(), [&] (Block& a) { splitterQueue.push_back(&a); });
+            
+            // Then perform the actual splitting until there are no more splitters.
+            while (!splitterQueue.empty()) {
+                refinePartition(backwardTransitions, *splitterQueue.front(), partition, splitterQueue);
+                splitterQueue.pop_front();
+            }
+            
+            // 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& block : partition.getBlocks()) {
+                // We need to sort the states to allow for rapid construction of the blocks.
+                std::sort(partition.getBeginOfStates(block), partition.getEndOfStates(block));
+                this->blocks[block.getId()] = block_type(partition.getBeginOfStates(block), partition.getEndOfStates(block), true);
+            }
+
+            // If we are required to build the quotient model, do so now.
+            if (buildQuotient) {
+                this->buildQuotient(model, partition);
+            }
+            
+            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+            
+            if (storm::settings::generalSettings().isShowStatisticsSet()) {
+                std::cout << "Computed bisimulation quotient in " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms." << std::endl;
+            }
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::refineBlockProbabilities(Block& block, Partition& partition, std::deque<Block*>& splitterQueue) {
+            // Sort the states in the block based on their probabilities.
+            std::sort(partition.getBeginOfStates(block), partition.getEndOfStates(block), [&partition] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return partition.getValue(a) < partition.getValue(b); } );
+            
+            // FIXME: This can probably be done more efficiently.
+            std::sort(partition.getBeginOfValues(block), partition.getEndOfValues(block));
+            
+            // Update the positions vector.
+            storm::storage::sparse::state_type position = block.getBegin();
+            for (auto stateIt = partition.getBeginOfStates(block), stateIte = partition.getEndOfStates(block); stateIt != stateIte; ++stateIt, ++position) {
+                partition.setPosition(*stateIt, position);
+            }
+            
+            // Finally, we need to scan the ranges of states that agree on the probability.
+            storm::storage::sparse::state_type beginIndex = block.getBegin();
+            storm::storage::sparse::state_type currentIndex = beginIndex;
+            storm::storage::sparse::state_type endIndex = block.getEnd();
+            
+            // 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(partition.getValueAtPosition(beginIndex), partition.getValueAtPosition(endIndex - 1))) {
+                // 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 = partition.getValueAtPosition(beginIndex);
+                
+                typename std::vector<ValueType>::const_iterator valueIterator = partition.getValues().begin() + beginIndex + 1;
+                ++currentIndex;
+                while (currentIndex < endIndex && comparator.isEqual(*valueIterator, currentValue)) {
+                    ++valueIterator;
+                    ++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();
+                }
+                beginIndex = currentIndex;
+            }
+        }
+        
+        template<typename ValueType>
+        void DeterministicModelStrongBisimulationDecomposition<ValueType>::refinePartition(storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Block& splitter, Partition& partition, std::deque<Block*>& splitterQueue) {
+            std::list<Block*> predecessorBlocks;
+            
+            // Iterate over all states of the splitter and check its predecessors.
+            storm::storage::sparse::state_type currentPosition = splitter.getBegin();
+            for (auto stateIterator = partition.getBeginOfStates(splitter), stateIte = partition.getEndOfStates(splitter); stateIterator != stateIte; ++stateIterator, ++currentPosition) {
+                storm::storage::sparse::state_type currentState = *stateIterator;
+                
+                uint_fast64_t elementsToSkip = 0;
+                for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) {
+                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+                    
+                    Block& predecessorBlock = partition.getBlock(predecessor);
+                    
+                    // If the predecessor block has just one state, there is no point in splitting 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.getStates().begin() + splitter.getOriginalBegin(), stateIte = partition.getStates().begin() + splitter.getMarkedPosition(); stateIterator != stateIte; ++stateIterator) {
+                storm::storage::sparse::state_type currentState = *stateIterator;
+                
+                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();
+                    }
+                }
+            }
+            
+            // Reset the marked position of the splitter to the begin.
+            splitter.setMarkedPosition(splitter.getBegin());
+            
+            std::list<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;
+                }
+                
+                refineBlockProbabilities(*blockPtr, partition, splitterQueue);
+            }
+        }
+        
+        template<typename ValueType>
+        template<typename ModelType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition DeterministicModelStrongBisimulationDecomposition<ValueType>::getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::string const& phiLabel, std::string const& psiLabel, bool bounded) {
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiLabel == "true" ? storm::storage::BitVector(model.getNumberOfStates(), true) : model.getLabeledStates(phiLabel), model.getLabeledStates(psiLabel));
+            Partition partition(model.getNumberOfStates(), statesWithProbability01.first, bounded ? model.getLabeledStates(psiLabel) : statesWithProbability01.second, phiLabel, psiLabel);
+            return partition;
+        }
+        
+        template<typename ValueType>
+        template<typename ModelType>
+        typename DeterministicModelStrongBisimulationDecomposition<ValueType>::Partition DeterministicModelStrongBisimulationDecomposition<ValueType>::getLabelBasedInitialPartition(ModelType const& model) {
+            Partition partition(model.getNumberOfStates());
+            for (auto const& label : model.getStateLabeling().getAtomicPropositions()) {
+                if (label == "init") {
+                    continue;
+                }
+                partition.splitLabel(model.getLabeledStates(label));
+            }
+            return partition;
+        }
+        
+        template class DeterministicModelStrongBisimulationDecomposition<double>;
+    }
+}
\ No newline at end of file
diff --git a/src/storage/DeterministicModelStrongBisimulationDecomposition.h b/src/storage/DeterministicModelStrongBisimulationDecomposition.h
new file mode 100644
index 000000000..e5ac547db
--- /dev/null
+++ b/src/storage/DeterministicModelStrongBisimulationDecomposition.h
@@ -0,0 +1,452 @@
+#ifndef STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_
+#define STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_
+
+#include <queue>
+#include <deque>
+
+#include "src/storage/sparse/StateType.h"
+#include "src/storage/Decomposition.h"
+#include "src/models/Dtmc.h"
+#include "src/models/Ctmc.h"
+#include "src/storage/Distribution.h"
+#include "src/utility/ConstantsComparator.h"
+#include "src/utility/OsDetection.h"
+
+namespace storm {
+    namespace storage {
+        
+        /*!
+         * This class represents the decomposition model into its (strong) bisimulation quotient.
+         */
+        template <typename ValueType>
+        class DeterministicModelStrongBisimulationDecomposition : public Decomposition<StateBlock> {
+        public:
+            /*!
+             * Decomposes the given DTMC into equivalence classes under strong bisimulation.
+             *
+             * @param model The model to decompose.
+             * @param buildQuotient Sets whether or not the quotient model is to be built.
+             */
+            DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc<ValueType> const& model, bool buildQuotient = false);
+
+            /*!
+             * Decomposes the given CTMC into equivalence classes under strong bisimulation.
+             *
+             * @param model The model to decompose.
+             * @param buildQuotient Sets whether or not the quotient model is to be built.
+             */
+            DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc<ValueType> const& model, bool buildQuotient = false);
+            
+            /*!
+             * Decomposes the given DTMC into equivalence classes under strong bisimulation in a way that onle safely
+             * preserves formulas of the form phi until psi.
+             *
+             * @param model The model to decompose.
+             * @param phiLabel The label that all phi states carry in the model.
+             * @param psiLabel The label that all psi states carry in the model.
+             * @param bounded If set to true, also bounded until formulas are preserved.
+             * @param buildQuotient Sets whether or not the quotient model is to be built.
+             */
+            DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc<ValueType> const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient = false);
+            
+            /*!
+             * Decomposes the given CTMC into equivalence classes under strong bisimulation in a way that onle safely
+             * preserves formulas of the form phi until psi.
+             *
+             * @param model The model to decompose.
+             * @param phiLabel The label that all phi states carry in the model.
+             * @param psiLabel The label that all psi states carry in the model.
+             * @param bounded If set to true, also bounded until formulas are preserved.
+             * @param buildQuotient Sets whether or not the quotient model is to be built.
+             */
+            DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc<ValueType> const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient = false);
+            
+            /*!
+             * Retrieves the quotient of the model under the previously computed bisimulation.
+             *
+             * @return The quotient model.
+             */
+            std::shared_ptr<storm::models::AbstractDeterministicModel<ValueType>> getQuotient() const;
+            
+        private:
+            class Partition;
+            
+            class Block {
+            public:
+                typedef typename std::list<Block>::const_iterator const_iterator;
+                
+                Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::shared_ptr<std::string> const& label = nullptr);
+                
+                // 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);
+
+                // Moves the end index of the block one step to the front.
+                void decrementEnd();
+                
+                // 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.
+                const_iterator getIterator() const;
+
+                // Returns the iterator the block in the list of overall blocks.
+                void setIterator(const_iterator it);
+                
+                // Returns the iterator the next block in the list of overall blocks if it exists.
+                const_iterator getNextIterator() const;
+
+                // Returns the iterator the next block in the list of overall blocks if it exists.
+                const_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;
+                
+                // Retrieves whether the block has a special label.
+                bool hasLabel() const;
+                
+                // Retrieves the special label of the block if it has one.
+                std::string const& getLabel() const;
+                
+                // Retrieves a pointer to the label of the block (which is the nullptr if there is none).
+                std::shared_ptr<std::string> const& getLabelPtr() const;
+                
+            private:
+                // An iterator to itself. This is needed to conveniently insert elements in the overall list of blocks
+                // kept by the partition.
+                const_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;
+                
+                // The label of the block or nullptr if it has none.
+                std::shared_ptr<std::string> label;
+                
+                // A counter for the IDs of the blocks.
+                static std::size_t blockId;
+            };
+            
+            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.
+                 */
+                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 otherLabel The label that is to be attached to all other states.
+                 * @param prob1Label The label that is to be attached to all states with probability 1.
+                 */
+                Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, std::string const& otherLabel, std::string const& prob1Label);
+                
+                Partition() = default;
+                Partition(Partition const& other) = default;
+                Partition& operator=(Partition const& other) = default;
+#ifndef WINDOWS
+                Partition(Partition&& other) = default;
+                Partition& operator=(Partition&& other) = default;
+#endif
+                
+                /*!
+                 * 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();
+                
+                // Retrieves the vector of all the states.
+                std::vector<storm::storage::sparse::state_type>& getStates();
+
+                // 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 getBeginOfStates(Block const& block);
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::iterator getEndOfStates(Block const& block);
+
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::const_iterator getBeginOfStates(Block const& block) const;
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                std::vector<storm::storage::sparse::state_type>::const_iterator getEndOfStates(Block const& block) const;
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                typename std::vector<ValueType>::iterator getBeginOfValues(Block const& block);
+                
+                // Returns an iterator to the beginning of the states of the given block.
+                typename std::vector<ValueType>::iterator getEndOfValues(Block const& block);
+
+                // 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<ValueType>& getValues();
+
+                // 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, std::vector<storm::storage::sparse::state_type>::iterator first, std::vector<storm::storage::sparse::state_type>::iterator end);
+                
+                // Retrieves the first block of the partition.
+                Block& getFirstBlock();
+            private:
+                // The list of blocks in the partition.
+                std::list<Block> blocks;
+                
+                // A mapping of states to their blocks.
+                std::vector<Block*> stateToBlockMapping;
+                
+                // A vector containing all the states. It is ordered in a special 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;
+                
+                // This vector stores the probabilities of going to the current splitter.
+                std::vector<ValueType> values;
+            };
+            
+            /*!
+             * 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 backwardTransitions The backward transitions of the model.
+             * @param The initial partition.
+             * @param buildQuotient If set, the quotient model is built and may be retrieved using the getQuotient()
+             * method.
+             */
+            template<typename ModelType>
+            void partitionRefinement(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Partition& partition, bool buildQuotient);
+            
+            /*!
+             * Refines the partition based on the provided splitter. After calling this method all blocks are stable
+             * with respect to the splitter.
+             *
+             * @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 splitterQueue A queue into which all blocks that were split are inserted so they can be treated
+             * as splitters in the future.
+             */
+            void refinePartition(storm::storage::SparseMatrix<ValueType> const& backwardTransitions, Block& splitter, Partition& partition, std::deque<Block*>& splitterQueue);
+            
+            /*!
+             * 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 splitterQueue A queue into which all blocks that were split are inserted so they can be treated
+             * as splitters in the future.
+             */
+            void refineBlockProbabilities(Block& block, Partition& partition, std::deque<Block*>& splitterQueue);
+            
+            /*!
+             * 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 partition The previously computed partition. This is used for quickly retrieving the block of a
+             * state.
+             */
+            template<typename ModelType>
+            void buildQuotient(ModelType const& model, Partition const& partition);
+
+            /*!
+             * 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 phiLabel The label that all phi states carry in the model.
+             * @param psiLabel The label that all psi states carry in the model.
+             * @param bounded If set to true, the initial partition will be chosen in such a way that preserves bounded
+             * reachability queries.
+             * @return The resulting partition.
+             */
+            template<typename ModelType>
+            Partition getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::string const& phiLabel, std::string const& psiLabel, bool bounded = false);
+            
+            /*!
+             * 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.
+             * @return The resulting partition.
+             */
+            template<typename ModelType>
+            Partition getLabelBasedInitialPartition(ModelType const& model);
+            
+            // If required, a quotient model is built and stored in this member.
+            std::shared_ptr<storm::models::AbstractDeterministicModel<ValueType>> quotient;
+            
+            // A comparator that is used for determining whether two probabilities are considered to be equal.
+            storm::utility::ConstantsComparator<ValueType> comparator;
+        };
+    }
+}
+
+#endif /* STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_ */
\ No newline at end of file
diff --git a/src/storage/Distribution.cpp b/src/storage/Distribution.cpp
new file mode 100644
index 000000000..50e93f567
--- /dev/null
+++ b/src/storage/Distribution.cpp
@@ -0,0 +1,100 @@
+#include "src/storage/Distribution.h"
+
+#include <algorithm>
+#include <iostream>
+
+#include "src/settings/SettingsManager.h"
+
+namespace storm {
+    namespace storage {
+        
+        template<typename ValueType>
+        Distribution<ValueType>::Distribution() : hash(0) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        bool Distribution<ValueType>::operator==(Distribution<ValueType> const& other) 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()) {
+                return false;
+            }
+            
+            auto first1 = this->distribution.begin();
+            auto last1 = this->distribution.end();
+            auto first2 = other.distribution.begin();
+            auto last2 = other.distribution.end();
+            
+            for (; first1 != last1; ++first1, ++first2) {
+                if (first1->first != first2->first) {
+                    return false;
+                }
+                if (std::abs(first1->second - first2->second) > 1e-6) {
+                    return false;
+                }
+            }
+            
+            return true;
+        }
+        
+        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);
+            }
+            this->distribution[state] += probability;
+        }
+        
+        template<typename ValueType>
+        typename Distribution<ValueType>::iterator Distribution<ValueType>::begin() {
+            return this->distribution.begin();
+        }
+
+        template<typename ValueType>
+        typename Distribution<ValueType>::const_iterator Distribution<ValueType>::begin() const {
+            return this->distribution.begin();
+        }
+        
+        template<typename ValueType>
+        typename Distribution<ValueType>::iterator Distribution<ValueType>::end() {
+            return this->distribution.end();
+        }
+        
+        template<typename ValueType>
+        typename Distribution<ValueType>::const_iterator Distribution<ValueType>::end() const {
+            return this->distribution.end();
+        }
+        
+        template<typename ValueType>
+        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;
+                this->distribution.erase(probabilityIterator);
+                
+                for (auto& entry : this->distribution) {
+                    entry.second *= scaleValue;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        std::size_t Distribution<ValueType>::getHash() const {
+            return this->hash ^ (this->distribution.size() << 8);
+        }
+        
+        template<typename ValueType>
+        std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution) {
+            out << "{";
+            for (auto const& entry : distribution) {
+                out << "[" << entry.second << ": " << entry.first << "], ";
+            }
+            out << "}";
+            
+            return out;
+        }
+        
+        template class Distribution<double>;
+        template std::ostream& operator<<(std::ostream& out, Distribution<double> const& distribution);
+    }
+}
diff --git a/src/storage/Distribution.h b/src/storage/Distribution.h
new file mode 100644
index 000000000..010cf43bc
--- /dev/null
+++ b/src/storage/Distribution.h
@@ -0,0 +1,110 @@
+#ifndef STORM_STORAGE_DISTRIBUTION_H_
+#define STORM_STORAGE_DISTRIBUTION_H_
+
+#include <vector>
+#include <ostream>
+#include <boost/container/flat_map.hpp>
+
+#include "src/storage/sparse/StateType.h"
+
+namespace storm {
+    namespace storage {
+        
+        template<typename ValueType>
+        class Distribution {
+        public:
+            typedef boost::container::flat_map<storm::storage::sparse::state_type, ValueType> container_type;
+            typedef typename container_type::iterator iterator;
+            typedef typename container_type::const_iterator const_iterator;
+            
+            /*!
+             * Creates an empty distribution.
+             */
+            Distribution();
+            
+            /*!
+             * Checks whether the two distributions specify the same probabilities to go to the same states.
+             *
+             * @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;
+            
+            /*!
+             * Assigns the given state the given probability under this distribution.
+             *
+             * @param state The state to which to assign the probability.
+             * @param probability The probability to assign.
+             */
+            void addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability);
+            
+            /*!
+             * 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.
+             *
+             * @return The iterator to the elements in this distribution.
+             */
+            const_iterator begin() const;
+            
+            /*!
+             * Retrieves an iterator past the elements in this distribution.
+             *
+             * @return The iterator past the elements in this distribution.
+             */
+            iterator end();
+
+            /*!
+             * Retrieves an iterator past the elements in this distribution.
+             *
+             * @return The iterator past the elements in this distribution.
+             */
+            const_iterator end() const;
+            
+            /*!
+             * Scales the distribution by multiplying all the probabilities with 1/p where p is the probability of moving
+             * to the given state and sets the probability of moving to the given state to zero. If the probability is
+             * already zero, this operation has no effect.
+             *
+             * @param state The state whose associated probability is used to scale the distribution.
+             */
+            void scale(storm::storage::sparse::state_type const& state);
+            
+            /*!
+             * Retrieves the hash value of the distribution.
+             *
+             * @return The hash value of the distribution.
+             */
+            std::size_t getHash() 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);
+    }
+}
+
+namespace std {
+    
+    template <typename ValueType>
+    struct hash<storm::storage::Distribution<ValueType>> {
+        std::size_t operator()(storm::storage::Distribution<ValueType> const& distribution) const {
+            return (distribution.getHash());
+        }
+    };
+    
+}
+
+
+#endif /* STORM_STORAGE_DISTRIBUTION_H_ */
\ No newline at end of file
diff --git a/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp b/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp
new file mode 100644
index 000000000..bf41414de
--- /dev/null
+++ b/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp
@@ -0,0 +1,262 @@
+#include "src/storage/NaiveDeterministicModelBisimulationDecomposition.h"
+
+#include <unordered_map>
+#include <chrono>
+
+namespace storm {
+    namespace storage {
+        
+        template<typename ValueType>
+        NaiveDeterministicModelBisimulationDecomposition<ValueType>::NaiveDeterministicModelBisimulationDecomposition(storm::models::Dtmc<ValueType> const& dtmc, bool weak) {
+            computeBisimulationEquivalenceClasses(dtmc, weak);
+        }
+        
+        template<typename ValueType>
+        void NaiveDeterministicModelBisimulationDecomposition<ValueType>::computeBisimulationEquivalenceClasses(storm::models::Dtmc<ValueType> const& dtmc, bool weak) {
+            std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
+            // We start by computing the initial partition. In particular, we also keep a mapping of states to their blocks.
+            std::vector<std::size_t> stateToBlockMapping(dtmc.getNumberOfStates());
+            storm::storage::BitVector labeledStates = dtmc.getLabeledStates("observe0Greater1");
+            this->blocks.emplace_back(labeledStates.begin(), labeledStates.end());
+            std::for_each(labeledStates.begin(), labeledStates.end(), [&] (storm::storage::sparse::state_type const& state) { stateToBlockMapping[state] = 0; } );
+            labeledStates.complement();
+            this->blocks.emplace_back(labeledStates.begin(), labeledStates.end());
+            std::for_each(labeledStates.begin(), labeledStates.end(), [&] (storm::storage::sparse::state_type const& state) { stateToBlockMapping[state] = 1; } );
+            
+            // Check whether any of the blocks is empty and remove it.
+            auto eraseIterator = std::remove_if(this->blocks.begin(), this->blocks.end(), [] (typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type const& a) { return a.empty(); });
+            this->blocks.erase(eraseIterator, this->blocks.end());
+        
+            // Create empty distributions for the two initial blocks.
+            std::vector<storm::storage::Distribution<ValueType>> distributions(2);
+            
+            // Retrieve the backward transitions to allow for better checking of states that need to be re-examined.
+            storm::storage::SparseMatrix<ValueType> const& backwardTransitions = dtmc.getBackwardTransitions();
+            
+            // Initially, both blocks are potential splitters. A splitter is marked as a pair in which the first entry
+            // is the ID of the parent block of the splitter and the second entry is the block ID of the splitter itself.
+            std::deque<std::size_t> refinementQueue;
+            storm::storage::BitVector blocksInRefinementQueue(this->size());
+            for (uint_fast64_t index = 0; index < this->blocks.size(); ++index) {
+                refinementQueue.push_back(index);
+            }
+            
+            // As long as there are blocks to refine, well, refine them.
+            uint_fast64_t iteration = 0;
+            while (!refinementQueue.empty()) {
+                ++iteration;
+                
+                // Optionally sort the queue to potentially speed up the convergence.
+                // std::sort(refinementQueue.begin(), refinementQueue.end(), [=] (std::size_t const& a, std::size_t const& b) { return this->blocks[a].size() > this->blocks[b].size(); });
+                
+                std::size_t currentBlock = refinementQueue.front();
+                refinementQueue.pop_front();
+                blocksInRefinementQueue.set(currentBlock, false);
+                
+                splitBlock(dtmc, backwardTransitions, currentBlock, stateToBlockMapping, distributions, blocksInRefinementQueue, refinementQueue, weak);
+            }
+                        
+            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+            
+            std::cout << "Bisimulation quotient has " << this->blocks.size() << " blocks and took " << iteration << " iterations and " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms." << std::endl;
+        }
+        
+        template<typename ValueType>
+        std::size_t NaiveDeterministicModelBisimulationDecomposition<ValueType>::splitPredecessorsGraphBased(storm::models::Dtmc<ValueType> const& dtmc, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::size_t const& blockId, std::vector<std::size_t>& stateToBlockMapping, std::vector<storm::storage::Distribution<ValueType>>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque<std::size_t>& graphRefinementQueue, storm::storage::BitVector& splitBlocks) {
+//            std::cout << "entering " << std::endl;
+            std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
+            
+//            std::cout << "refining predecessors of " << blockId << std::endl;
+            
+            // This method tries to split the blocks of predecessors of the provided block by checking whether there is
+            // a transition into the current block or not.
+            std::unordered_map<std::size_t, typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type> predecessorBlockToNewBlock;
+            
+            // Now for each predecessor block which state could actually reach the current block.
+            std::chrono::high_resolution_clock::time_point predIterStart = std::chrono::high_resolution_clock::now();
+            int predCounter = 0;
+            storm::storage::BitVector preds(dtmc.getNumberOfStates());
+            for (auto const& state : this->blocks[blockId]) {
+                for (auto const& predecessorEntry : backwardTransitions.getRow(state)) {
+                    ++predCounter;
+                    storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+                    if (!preds.get(predecessor)) {
+//                        std::cout << "having pred " << predecessor << " with block " << stateToBlockMapping[predecessor] << std::endl;
+                        predecessorBlockToNewBlock[stateToBlockMapping[predecessor]].insert(predecessor);
+                        preds.set(predecessor);
+                    }
+                }
+            }
+            std::chrono::high_resolution_clock::duration predIterTime = std::chrono::high_resolution_clock::now() - predIterStart;
+            std::cout << "took " << std::chrono::duration_cast<std::chrono::milliseconds>(predIterTime).count() << "ms. for " << predCounter << " predecessors" << std::endl;
+
+            // Now, we can check for each predecessor block whether it needs to be split.
+            for (auto const& blockNewBlockPair : predecessorBlockToNewBlock) {
+                if (this->blocks[blockNewBlockPair.first].size() > blockNewBlockPair.second.size()) {
+                    std::cout << "splitting!" << std::endl;
+//                    std::cout << "found that not all " << this->blocks[blockNewBlockPair.first].size() << " states of block " << blockNewBlockPair.first << " have a successor in " << blockId << " but only " << blockNewBlockPair.second.size() << std::endl;
+//                    std::cout << "original: " << this->blocks[blockNewBlockPair.first] << std::endl;
+//                    std::cout << "states with pred: " << blockNewBlockPair.second << std::endl;
+                    // Now update the block mapping for the smaller of the two blocks.
+                    typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type smallerBlock;
+                    typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type largerBlock;
+                    if (blockNewBlockPair.second.size() < this->blocks[blockNewBlockPair.first].size()/2) {
+                        smallerBlock = std::move(blockNewBlockPair.second);
+                        std::set_difference(this->blocks[blockNewBlockPair.first].begin(), this->blocks[blockNewBlockPair.first].end(), smallerBlock.begin(), smallerBlock.end(), std::inserter(largerBlock, largerBlock.begin()));
+                    } else {
+                        largerBlock = std::move(blockNewBlockPair.second);
+                        std::set_difference(this->blocks[blockNewBlockPair.first].begin(), this->blocks[blockNewBlockPair.first].end(), largerBlock.begin(), largerBlock.end(), std::inserter(smallerBlock, smallerBlock.begin()));
+                    }
+                    
+//                    std::cout << "created a smaller block with " << smallerBlock.size() << " states and a larger one with " << largerBlock.size() << "states" << std::endl;
+//                    std::cout << "smaller: " << smallerBlock << std::endl;
+//                    std::cout << "larger: " << largerBlock << std::endl;
+                    
+                    this->blocks.emplace_back(std::move(smallerBlock));
+                    this->blocks[blockNewBlockPair.first] = std::move(largerBlock);
+                    
+                    // Update the block mapping of all moved states.
+                    std::size_t newBlockId = this->blocks.size() - 1;
+                    for (auto const& state : this->blocks.back()) {
+                        stateToBlockMapping[state] = newBlockId;
+//                        std::cout << "updating " << state << " to block " << newBlockId << std::endl;
+                    }
+                    
+                    blocksInRefinementQueue.resize(this->blocks.size());
+                    
+                    // Add both parts of the old block to the queue.
+                    if (!blocksInRefinementQueue.get(blockNewBlockPair.first)) {
+                        graphRefinementQueue.push_back(blockNewBlockPair.first);
+                        blocksInRefinementQueue.set(blockNewBlockPair.first, true);
+//                        std::cout << "adding " << blockNewBlockPair.first << " to refine further using graph-based analysis " << std::endl;
+                    }
+                    
+                    graphRefinementQueue.push_back(this->blocks.size() - 1);
+                    blocksInRefinementQueue.set(this->blocks.size() - 1);
+//                    std::cout << "adding " << this->blocks.size() - 1 << " to refine further using graph-based analysis " << std::endl;
+                }
+            }
+            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+            std::cout << "refinement of predecessors of block " << blockId << " took " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms." << std::endl;
+
+//            std::cout << "done. "<< std::endl;
+            
+            // FIXME
+            return 0;
+        }
+        
+        template<typename ValueType>
+        std::size_t NaiveDeterministicModelBisimulationDecomposition<ValueType>::splitBlock(storm::models::Dtmc<ValueType> const& dtmc, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::size_t const& blockId, std::vector<std::size_t>& stateToBlockMapping, std::vector<storm::storage::Distribution<ValueType>>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque<std::size_t>& refinementQueue, bool weakBisimulation) {
+            std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
+            std::unordered_map<storm::storage::Distribution<ValueType>, typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type> distributionToNewBlocks;
+            
+            // Traverse all states of the block and check whether they have different distributions.
+            std::chrono::high_resolution_clock::time_point gatherStart = std::chrono::high_resolution_clock::now();
+            for (auto const& state : this->blocks[blockId]) {
+                // Now construct the distribution of this state wrt. to the current partitioning.
+                storm::storage::Distribution<ValueType> distribution;
+                for (auto const& successorEntry : dtmc.getTransitionMatrix().getRow(state)) {
+                    distribution.addProbability(static_cast<storm::storage::sparse::state_type>(stateToBlockMapping[successorEntry.getColumn()]), successorEntry.getValue());
+                }
+                
+                // If we are requested to compute a weak bisimulation, we need to scale the distribution with the
+                // self-loop probability.
+                if (weakBisimulation) {
+                    distribution.scale(blockId);
+                }
+                
+                // If the distribution already exists, we simply add the state. Otherwise, we open a new block.
+                auto distributionIterator = distributionToNewBlocks.find(distribution);
+                if (distributionIterator != distributionToNewBlocks.end()) {
+                    distributionIterator->second.insert(state);
+                } else {
+                    distributionToNewBlocks[distribution].insert(state);
+                }
+            }
+            
+            std::chrono::high_resolution_clock::duration gatherTime = std::chrono::high_resolution_clock::now() - gatherStart;
+            std::cout << "time to iterate over all states was " << std::chrono::duration_cast<std::chrono::milliseconds>(gatherTime).count() << "ms." << std::endl;
+            
+            // Now we are ready to split the block.
+            if (distributionToNewBlocks.size() == 1) {
+                // If there is just one behavior, we just set the distribution as the new one for this block.
+                // distributions[blockId] = std::move(distributionToNewBlocks.begin()->first);
+            } else {
+                // In this case, we need to split the block.
+                typename NaiveDeterministicModelBisimulationDecomposition<ValueType>::block_type tmpBlock;
+                
+                auto distributionIterator = distributionToNewBlocks.begin();
+                STORM_LOG_ASSERT(distributionIterator != distributionToNewBlocks.end(), "One block has no distribution...");
+
+                //                distributions[blockId] = std::move(distributionIterator->first);
+                tmpBlock = std::move(distributionIterator->second);
+                std::swap(this->blocks[blockId], tmpBlock);
+                ++distributionIterator;
+                
+                // Remember the number of blocks prior to splitting for later use.
+                std::size_t beforeNumberOfBlocks = this->blocks.size();
+                
+                for (; distributionIterator != distributionToNewBlocks.end(); ++distributionIterator) {
+                    // In this case, we need to move the newly created block to the end of the list of actual blocks.
+                    this->blocks.emplace_back(std::move(distributionIterator->second));
+                    distributions.emplace_back(std::move(distributionIterator->first));
+                    
+                    // Update the mapping of states to their blocks.
+                    std::size_t newBlockId = this->blocks.size() - 1;
+                    for (auto const& state : this->blocks.back()) {
+                        stateToBlockMapping[state] = newBlockId;
+                    }
+                }
+                
+                // Now that we have split the block, we try to trigger a chain of graph-based splittings. That is, we
+                // try to split the predecessors of the current block by checking whether they have some transition
+                // to one given sub-block of the current block.
+                std::deque<std::size_t> localRefinementQueue;
+                storm::storage::BitVector blocksInLocalRefinementQueue(this->size());
+                localRefinementQueue.push_back(blockId);
+//                std::cout << "adding " << blockId << " to local ref queue " << std::endl;
+                blocksInLocalRefinementQueue.set(blockId);
+                for (std::size_t i = beforeNumberOfBlocks; i < this->blocks.size(); ++i) {
+                    localRefinementQueue.push_back(i);
+                    blocksInLocalRefinementQueue.set(i);
+//                    std::cout << "adding " << i << " to local refinement queue " << std::endl;
+                }
+
+                // We need to keep track of which blocks were split so we can later add all their predecessors as
+                // candidates for furter splitting.
+                storm::storage::BitVector splitBlocks(this->size());
+                
+//                while (!localRefinementQueue.empty()) {
+//                    std::size_t currentBlock = localRefinementQueue.front();
+//                    localRefinementQueue.pop_front();
+//                    blocksInLocalRefinementQueue.set(currentBlock, false);
+//                    
+//                    splitPredecessorsGraphBased(dtmc, backwardTransitions, currentBlock, stateToBlockMapping, distributions, blocksInLocalRefinementQueue, localRefinementQueue, splitBlocks);
+//                }
+                
+//                std::cout << "split blocks: " << std::endl;
+//                std::cout << splitBlocks << std::endl;
+                
+                // Since we created some new blocks, we need to extend the bit vector storing the blocks in the
+                // refinement queue.
+                blocksInRefinementQueue.resize(this->blocks.size());
+                
+                // Insert blocks that possibly need a refinement into the queue.
+                for (auto const& state : tmpBlock) {
+                    for (auto const& predecessor : backwardTransitions.getRow(state)) {
+                        if (!blocksInRefinementQueue.get(stateToBlockMapping[predecessor.getColumn()])) {
+                            blocksInRefinementQueue.set(stateToBlockMapping[predecessor.getColumn()]);
+                            refinementQueue.push_back(stateToBlockMapping[predecessor.getColumn()]);
+                        }
+                    }
+                }
+            }
+            
+            std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+            std::cout << "refinement of block " << blockId << " took " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms." << std::endl;
+            return distributionToNewBlocks.size();
+        }
+        
+        template class NaiveDeterministicModelBisimulationDecomposition<double>;
+    }
+}
\ No newline at end of file
diff --git a/src/storage/NaiveDeterministicModelBisimulationDecomposition.h b/src/storage/NaiveDeterministicModelBisimulationDecomposition.h
new file mode 100644
index 000000000..300922cf5
--- /dev/null
+++ b/src/storage/NaiveDeterministicModelBisimulationDecomposition.h
@@ -0,0 +1,35 @@
+#ifndef STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+#define STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_
+
+#include <queue>
+#include <deque>
+
+#include "src/storage/Decomposition.h"
+#include "src/models/Dtmc.h"
+#include "src/storage/Distribution.h"
+
+namespace storm {
+    namespace storage {
+        
+        /*!
+         * This class represents the decomposition model into its bisimulation quotient.
+         */
+        template <typename ValueType>
+        class NaiveDeterministicModelBisimulationDecomposition : public Decomposition<StateBlock> {
+        public:
+            NaiveDeterministicModelBisimulationDecomposition() = default;
+            
+            /*!
+             * Decomposes the given DTMC into equivalence classes under weak or strong bisimulation.
+             */
+            NaiveDeterministicModelBisimulationDecomposition(storm::models::Dtmc<ValueType> const& model, bool weak = false);
+            
+        private:
+            void computeBisimulationEquivalenceClasses(storm::models::Dtmc<ValueType> const& model, bool weak);
+            std::size_t splitPredecessorsGraphBased(storm::models::Dtmc<ValueType> const& dtmc, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::size_t const& block, std::vector<std::size_t>& stateToBlockMapping, std::vector<storm::storage::Distribution<ValueType>>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque<std::size_t>& refinementQueue, storm::storage::BitVector& splitBlocks);
+            std::size_t splitBlock(storm::models::Dtmc<ValueType> const& dtmc, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::size_t const& block, std::vector<std::size_t>& stateToBlockMapping, std::vector<storm::storage::Distribution<ValueType>>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque<std::size_t>& refinementQueue, bool weakBisimulation);
+        };
+    }
+}
+
+#endif /* STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ */
\ No newline at end of file
diff --git a/src/storage/StateBlock.cpp b/src/storage/StateBlock.cpp
index fb4f21133..bc270872d 100644
--- a/src/storage/StateBlock.cpp
+++ b/src/storage/StateBlock.cpp
@@ -30,6 +30,10 @@ namespace storm {
             states.insert(state);
         }
 
+        StateBlock::iterator StateBlock::insert(container_type::const_iterator iterator, value_type const& state) {
+            return states.insert(iterator, state);
+        }
+        
         void StateBlock::erase(value_type const& state) {
             states.erase(state);
         }
@@ -55,5 +59,6 @@ namespace storm {
             out << block.getStates();
             return out;
         }
+        
     }
 }
\ No newline at end of file
diff --git a/src/storage/StateBlock.h b/src/storage/StateBlock.h
index dcd516277..e9a8d5ebe 100644
--- a/src/storage/StateBlock.h
+++ b/src/storage/StateBlock.h
@@ -2,6 +2,7 @@
 #define STORM_STORAGE_BLOCK_H_
 
 #include <boost/container/flat_set.hpp>
+#include <boost/container/container_fwd.hpp>
 
 #include "src/utility/OsDetection.h"
 #include "src/storage/sparse/StateType.h"
@@ -35,10 +36,15 @@ namespace storm {
              *
              * @param first The first element of the range to insert.
              * @param last The last element of the range (that is itself not inserted).
+             * @param sortedAndUnique If set to true, the input range is assumed to be sorted and duplicate-free.
              */
             template <typename InputIterator>
-            StateBlock(InputIterator first, InputIterator last) : states(first, last) {
-                // Intentionally left empty.
+            StateBlock(InputIterator first, InputIterator last, bool sortedAndUnique = false) {
+                if (sortedAndUnique) {
+                    this->states = container_type(boost::container::ordered_unique_range_t(), first, last);
+                } else {
+                    this->states = container_type(first, last);
+                }
             }
             
             /*!
@@ -101,6 +107,13 @@ namespace storm {
              * @param state The state to add to this SCC.
              */
             void insert(value_type const& state);
+
+            /*!
+             * Inserts the given element into this SCC.
+             *
+             * @param state The state to add to this SCC.
+             */
+            iterator insert(container_type::const_iterator iterator, value_type const& state);
             
             /*!
              * Removes the given element from this SCC.
diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h
new file mode 100644
index 000000000..fc46b251f
--- /dev/null
+++ b/src/utility/ConstantsComparator.h
@@ -0,0 +1,84 @@
+#ifndef STORM_UTILITY_CONSTANTSCOMPARATOR_H_
+#define STORM_UTILITY_CONSTANTSCOMPARATOR_H_
+
+#ifdef max
+#	undef max
+#endif
+
+#ifdef min
+#	undef min
+#endif
+
+#include <limits>
+#include <cstdint>
+
+#include "src/settings/SettingsManager.h"
+
+namespace storm {
+    namespace utility {
+        
+        template<typename ValueType>
+        ValueType one() {
+            return ValueType(1);
+        }
+        
+        template<typename ValueType>
+        ValueType zero() {
+            return ValueType(0);
+        }
+        
+        template<typename ValueType>
+        ValueType infinity() {
+            return std::numeric_limits<ValueType>::infinity();
+        }
+        
+        // A class that can be used for comparing constants.
+        template<typename ValueType>
+        class ConstantsComparator {
+        public:
+            bool isOne(ValueType const& value) {
+                return value == one<ValueType>();
+            }
+            
+            bool isZero(ValueType const& value) {
+                return value == zero<ValueType>();
+            }
+            
+            bool isEqual(ValueType const& value1, ValueType const& value2) {
+                return value1 == value2;
+            }
+        };
+        
+        // For doubles we specialize this class and consider the comparison modulo some predefined precision.
+        template<>
+        class ConstantsComparator<double> {
+        public:
+            ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) {
+                // Intentionally left empty.
+            }
+
+            ConstantsComparator(double precision) : precision(precision) {
+                // Intentionally left empty.
+            }
+            
+            bool isOne(double const& value) {
+                return std::abs(value - one<double>()) <= precision;
+            }
+            
+            bool isZero(double const& value) {
+                return std::abs(value) <= precision;
+            }
+            
+            bool isEqual(double const& value1, double const& value2) {
+                return std::abs(value1 - value2) <= precision;
+            }
+            
+        private:
+            // The precision used for comparisons.
+            double precision;
+        };
+        
+    }
+}
+
+#endif /* STORM_UTILITY_CONSTANTSCOMPARATOR_H_ */
\ No newline at end of file
diff --git a/src/utility/cli.h b/src/utility/cli.h
index fd1c2a496..227e09dfb 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -8,6 +8,7 @@
 #include <sstream>
 #include <memory>
 
+#include "storm-config.h"
 // Includes for the linked libraries and versions header.
 #ifdef STORM_HAVE_INTELTBB
 #	include "tbb/tbb_stddef.h"
@@ -45,6 +46,10 @@ log4cplus::Logger printer;
 // Headers of adapters.
 #include "src/adapters/ExplicitModelAdapter.h"
 
+// Headers for model processing.
+#include "src/storage/NaiveDeterministicModelBisimulationDecomposition.h"
+#include "src/storage/DeterministicModelStrongBisimulationDecomposition.h"
+
 // Headers for counterexample generation.
 #include "src/counterexamples/MILPMinimalLabelSetGenerator.h"
 #include "src/counterexamples/SMTMinimalCommandSetGenerator.h"
@@ -258,6 +263,22 @@ namespace storm {
                     STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "No input model.");
                 }
                 
+                // Print some information about the model.
+                result->printModelInformationToStream(std::cout);
+                
+                if (settings.isBisimulationSet()) {
+                    STORM_LOG_THROW(result->getType() == storm::models::DTMC, storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only compatible with DTMCs.");
+                    std::shared_ptr<storm::models::Dtmc<double>> dtmc = result->template as<storm::models::Dtmc<double>>();
+                    
+                    STORM_PRINT(std::endl << "Performing bisimulation minimization..." << std::endl);
+                    storm::storage::DeterministicModelStrongBisimulationDecomposition<double> bisimulationDecomposition(*dtmc, true);
+                    
+                    result = bisimulationDecomposition.getQuotient();
+                    
+                    STORM_PRINT_AND_LOG(std::endl << "Model after minimization:" << std::endl);
+                    result->printModelInformationToStream(std::cout);
+                }
+                
                 return result;
             }
             
@@ -268,7 +289,7 @@ namespace storm {
                     
                     std::shared_ptr<storm::models::Mdp<double>> mdp = model->as<storm::models::Mdp<double>>();
 
-                    // FIXME: re-parse the program.
+                    // FIXME: do not re-parse the program.
                     std::string const programFile = storm::settings::generalSettings().getSymbolicModelFilename();
                     std::string const constants = storm::settings::generalSettings().getConstantDefinitionString();
                     storm::prism::Program program = storm::parser::PrismParser::parse(programFile);
@@ -297,15 +318,13 @@ namespace storm {
                 // Start by parsing/building the model.
                 std::shared_ptr<storm::models::AbstractModel<double>> model = buildModel<double>();
                 
-                // Print some information about the model.
-                model->printModelInformationToStream(std::cout);
-                
                 // If we were requested to generate a counterexample, we now do so.
                 if (settings.isCounterexampleSet()) {
                     STORM_LOG_THROW(settings.isPctlPropertySet(), storm::exceptions::InvalidSettingsException, "Unable to generate counterexample without a property.");
                     std::shared_ptr<storm::properties::prctl::PrctlFilter<double>> filterFormula = storm::parser::PrctlParser::parsePrctlFormula(settings.getPctlProperty());
                     generateCounterexample(model, filterFormula->getChild());
                 }
+                
             }
 
         }
diff --git a/src/utility/constants.h b/src/utility/constants.h
index 699ee7368..60e601aaa 100644
--- a/src/utility/constants.h
+++ b/src/utility/constants.h
@@ -254,8 +254,6 @@ inline bool isZero(storm::Polynomial const& p)
 }
 #endif
 
-
-
 template<typename T>
 inline bool isOne(T const& sum)
 {
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 33f637877..cbfad8758 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -204,6 +204,25 @@ namespace storm {
                 return result;
             }
             
+            /*!
+             * Computes the sets of states that have probability 0 or 1, respectively, of satisfying phi until psi in a
+             * deterministic model.
+             *
+             * @param backwardTransitions The backward transitions of the model whose graph structure to search.
+             * @param phiStates The set of all states satisfying phi.
+             * @param psiStates The set of all states satisfying psi.
+             * @return A pair of bit vectors such that the first bit vector stores the indices of all states
+             * with probability 0 and the second stores all indices of states with probability 1.
+             */
+            template <typename T>
+            static std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01(storm::storage::SparseMatrix<T> backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
+                result.first = performProbGreater0(backwardTransitions, phiStates, psiStates);
+                result.second = performProb1(backwardTransitions, phiStates, psiStates, result.first);
+                result.first.complement();
+                return result;
+            }
+            
             /*!
              * Computes the sets of states that have probability greater 0 of satisfying phi until psi under at least
              * one possible resolution of non-determinism in a non-deterministic model. Stated differently,