more work on reimplementation of weak bisim

Former-commit-id: 5bdd8ea139

src/storage/bisimulation/BisimulationDecomposition.cpp

@@ -160,6 +160,9 @@ namespace storm {
             this->performPartitionRefinement();
             std::chrono::high_resolution_clock::duration refinementTime = std::chrono::high_resolution_clock::now() - refinementStart;
             
+            std::cout << "final partition: " << std::endl;
+            this->partition.print();
+            
             std::chrono::high_resolution_clock::time_point extractionStart = std::chrono::high_resolution_clock::now();
             this->extractDecompositionBlocks();
             std::chrono::high_resolution_clock::duration extractionTime = std::chrono::high_resolution_clock::now() - extractionStart;

src/storage/bisimulation/Block.h

@@ -20,7 +20,7 @@ namespace storm {
                 
                 // Creates a new block with the given begin and end.
                 Block(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, Block* previous, Block* next, std::size_t id);
-                
+                                
                 Block() = default;
                 Block(Block const& other) = default;
                 Block& operator=(Block const& other) = default;

src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp

@@ -4,7 +4,7 @@
 #include <unordered_map>
 #include <chrono>
 #include <iomanip>
-#include <boost/iterator/transform_iterator.hpp>
+#include <boost/iterator/zip_iterator.hpp>
 
 #include "src/adapters/CarlAdapter.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
@@ -45,10 +45,12 @@ namespace storm {
             stateStack.reserve(this->model.getNumberOfStates());
             storm::storage::BitVector nondivergentStates(this->model.getNumberOfStates());
             
-            for (auto const& blockPtr : this->partition.getBlocks()) {
+            uint_fast64_t currentSize = this->partition.size();
+            for (uint_fast64_t blockIndex = 0; blockIndex < currentSize; ++blockIndex) {
+                auto& block = *this->partition.getBlocks()[blockIndex];
                 nondivergentStates.clear();
                 
-                for (auto stateIt = this->partition.begin(*blockPtr), stateIte = this->partition.end(*blockPtr); stateIt != stateIte; ++stateIt) {
+                for (auto stateIt = this->partition.begin(block), stateIte = this->partition.end(block); stateIt != stateIte; ++stateIt) {
                     if (nondivergentStates.get(*stateIt)) {
                         continue;
                     }
@@ -59,7 +61,7 @@ namespace storm {
                     for (auto const& successor : this->model.getRows(*stateIt)) {
                         // If there is such a transition, then we can mark all states in the current block that can
                         // reach the state as non-divergent.
-                        if (this->partition.getBlock(successor.getColumn()) != *blockPtr) {
+                        if (this->partition.getBlock(successor.getColumn()) != block) {
                             isDirectlyNonDivergent = true;
                             break;
                         }
@@ -74,7 +76,7 @@ namespace storm {
                             nondivergentStates.set(currentState);
                             
                             for (auto const& predecessor : this->backwardTransitions.getRow(currentState)) {
-                                if (this->partition.getBlock(predecessor.getColumn()) == *blockPtr && !nondivergentStates.get(predecessor.getColumn())) {
+                                if (this->partition.getBlock(predecessor.getColumn()) == block && !nondivergentStates.get(predecessor.getColumn())) {
                                     stateStack.push_back(predecessor.getColumn());
                                 }
                             }
@@ -82,18 +84,17 @@ namespace storm {
                     }
                 }
                 
-                
-                if (nondivergentStates.getNumberOfSetBits() > 0 && nondivergentStates.getNumberOfSetBits() != blockPtr->getNumberOfStates()) {
+                if (!nondivergentStates.empty() && nondivergentStates.getNumberOfSetBits() != block.getNumberOfStates()) {
                     // After performing the split, the current block will contain the divergent states only.
-                    this->partition.splitStates(*blockPtr, nondivergentStates);
+                    this->partition.splitStates(block, nondivergentStates);
                     
                     // Since the remaining states in the block are divergent, we can mark the block as absorbing.
                     // This also guarantees that the self-loop will be added to the state of the quotient
                     // representing this block of states.
-                    blockPtr->data().setAbsorbing(true);
-                } else if (nondivergentStates.getNumberOfSetBits() == 0) {
+                    block.data().setAbsorbing(true);
+                } else if (nondivergentStates.empty()) {
                     // If there are only diverging states in the block, we need to make it absorbing.
-                    blockPtr->data().setAbsorbing(true);
+                    block.data().setAbsorbing(true);
                 }
             }
         }
@@ -147,13 +148,18 @@ namespace storm {
             return this->comparator.isOne(silentProbabilities[state]);
         }
         
+        template<typename ModelType>
+        bool DeterministicModelBisimulationDecomposition<ModelType>::hasNonZeroSilentProbability(storm::storage::sparse::state_type const& state) const {
+            return !this->comparator.isZero(silentProbabilities[state]);
+        }
+        
         template<typename ModelType>
         typename DeterministicModelBisimulationDecomposition<ModelType>::ValueType DeterministicModelBisimulationDecomposition<ModelType>::getSilentProbability(storm::storage::sparse::state_type const& state) const {
             return silentProbabilities[state];
         }
         
         template<typename ModelType>
-        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitter(std::list<Block<BlockDataType>*>& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitterStrong(std::list<Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
             for (auto block : predecessorBlocks) {
                 // Depending on the actions we need to take, the block to refine changes, so we need to keep track of it.
                 Block<BlockDataType>* blockToRefineProbabilistically = block;
@@ -161,9 +167,8 @@ namespace storm {
                 bool split = false;
                 // If the new begin index has shifted to a non-trivial position, we need to split the block.
                 if (block->getBeginIndex() != block->data().marker1() && block->getEndIndex() != block->data().marker1()) {
-                    auto result = this->partition.splitBlock(*block, block->data().marker1());
-                    STORM_LOG_ASSERT(result.second, "Expected to split block, but that did not happen.");
                     split = true;
+                    this->partition.splitBlock(*block, block->data().marker1());
                     blockToRefineProbabilistically = block->getPreviousBlockPointer();
                 }
                 
@@ -243,7 +248,7 @@ namespace storm {
                     this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1);
                 }
                 
-                // Since we had to move an already explored state to the right of the current position, 
+                // Since we had to move an already explored state to the right of the current position,
                 ++elementsToSkip;
                 predecessorBlock.data().incrementMarker1();
                 predecessorBlock.data().incrementMarker2();
@@ -266,7 +271,7 @@ namespace storm {
                     
                     // We keep track of the probability of the predecessor moving to the splitter.
                     increaseProbabilityToSplitter(predecessor, predecessorBlock, predecessorEntry.getValue());
-
+                    
                     // Only move the state if it has not been seen as a predecessor before.
                     storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
                     if (predecessorPosition >= predecessorBlock.data().marker1()) {
@@ -279,6 +284,169 @@ namespace storm {
             splitter.data().setMarker2(splitter.getBeginIndex());
         }
         
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(bisimulation::Block<BlockDataType>& block) {
+            // For all states that do not have a successor in the block itself, we need to set the silent probability to 0.
+            std::for_each(silentProbabilities.begin() + block.data().marker1(), silentProbabilities.begin() + block.getEndIndex(), [] (ValueType& val) { val = storm::utility::zero<ValueType>(); } );
+            
+            // For all states that do have a successor in the block, we set the silent probability to the probability
+            // stored in the vector of probabilities going to the splitter.
+            auto it = boost::make_zip_iterator(boost::make_tuple(silentProbabilities.begin() + block.getBeginIndex(), probabilitiesToCurrentSplitter.begin() + block.getBeginIndex()));
+            auto ite = boost::make_zip_iterator(boost::make_tuple(silentProbabilities.begin() + block.data().marker1(), probabilitiesToCurrentSplitter.begin() + block.data().marker1()));
+            std::for_each(it, ite, [] (boost::tuple<ValueType&, ValueType&> const& tuple) { boost::get<0>(tuple) = boost::get<1>(tuple); } );
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::updateSilentProbabilitiesBasedOnTransitions(bisimulation::Block<BlockDataType>& block) {
+            for (auto stateIt = this->partition.begin(block), stateIte = this->partition.end(block); stateIt != stateIte; ++stateIt) {
+                if (hasNonZeroSilentProbability(*stateIt)) {
+                    ValueType newSilentProbability = storm::utility::zero<ValueType>();
+                    for (auto const& successorEntry : this->model.getTransitionMatrix().getRow(*stateIt)) {
+                        if (this->partition.getBlock(successorEntry.getColumn()) == block) {
+                            newSilentProbability += successorEntry.getValue();
+                        }
+                    }
+                    silentProbabilities[*stateIt] = newSilentProbability;
+                }
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::computeConditionalProbabilitiesForNonSilentStates(bisimulation::Block<BlockDataType>& block) {
+            auto it = boost::make_zip_iterator(boost::make_tuple(silentProbabilities.begin() + block.getBeginIndex(), probabilitiesToCurrentSplitter.begin() + block.getBeginIndex()));
+            auto ite = boost::make_zip_iterator(boost::make_tuple(silentProbabilities.begin() + block.data().marker1(), probabilitiesToCurrentSplitter.begin() + block.data().marker1()));
+            std::for_each(it, ite, [this] (boost::tuple<ValueType&, ValueType&> const& tuple) {
+                if (!this->comparator.isZero(boost::get<0>(tuple))) {
+                    boost::get<1>(tuple) /= storm::utility::one<ValueType>() - boost::get<0>(tuple);
+                }
+            } );
+        }
+        
+        template<typename ModelType>
+        std::vector<uint_fast64_t> DeterministicModelBisimulationDecomposition<ModelType>::computeNonSilentBlocks(bisimulation::Block<BlockDataType>& block) {
+            auto less = [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { return probabilitiesToCurrentSplitter[state1] < probabilitiesToCurrentSplitter[state2]; };
+            this->partition.sortRange(block.getBeginIndex(), block.data().marker1(), less);
+            return this->partition.computeRangesOfEqualValue(block.getBeginIndex(), block.data().marker1(), less);
+        }
+        
+        template<typename ModelType>
+        std::vector<storm::storage::BitVector> DeterministicModelBisimulationDecomposition<ModelType>::computeWeakStateLabelingBasedOnNonSilentBlocks(bisimulation::Block<BlockDataType> const& block, std::vector<uint_fast64_t> const& nonSilentBlockIndices) {
+            // Now that we have the split points of the non-silent states, we perform a backward search from
+            // each non-silent state and label the predecessors with the class of the non-silent state.
+            std::vector<storm::storage::BitVector> stateLabels(block.getNumberOfStates(), storm::storage::BitVector(nonSilentBlockIndices.size() - 1));
+            std::cout << "creating " << block.getNumberOfStates() << " labels " << std::endl;
+            
+            std::vector<storm::storage::sparse::state_type> stateStack;
+            stateStack.reserve(block.getNumberOfStates());
+            for (uint_fast64_t stateClassIndex = 0; stateClassIndex < nonSilentBlockIndices.size() - 1; ++stateClassIndex) {
+                for (auto stateIt = this->partition.begin() + nonSilentBlockIndices[stateClassIndex], stateIte = this->partition.begin() + nonSilentBlockIndices[stateClassIndex + 1]; stateIt != stateIte; ++stateIt) {
+                    std::cout << "moving backward from state " << *stateIt << std::endl;
+                    
+                    stateStack.push_back(*stateIt);
+                    stateLabels[this->partition.getPosition(*stateIt) - block.getBeginIndex()].set(stateClassIndex);
+                    while (!stateStack.empty()) {
+                        storm::storage::sparse::state_type currentState = stateStack.back();
+                        stateStack.pop_back();
+                        
+                        for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
+                            storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn();
+
+                            std::cout << " found pred " << predecessor << " with value " << predecessorEntry.getValue() << std::endl;
+                            if (this->comparator.isZero(predecessorEntry.getValue())) {
+                                continue;
+                            }
+                            
+                            // Only if the state is in the same block, is a silent state and it has not yet been
+                            // labeled with the current label.
+                            std::cout << "pred is at " << this->partition.getPosition(predecessor) << std::endl;
+                            std::cout << "begin of block " << block.getBeginIndex() << std::endl;
+                            std::cout << "accessing index " << (this->partition.getPosition(predecessor) - block.getBeginIndex()) << " of " << stateLabels.size() << std::endl;
+                            std::cout << predecessor << " is in the same block as " << block.getId() << " ? " << (this->partition.getBlock(predecessor) == block) << std::endl;
+//                            std::cout << isSilent(predecessor) << ", " << stateLabels[this->partition.getPosition(predecessor) - block.getBeginIndex()].get(stateClassIndex) << std::endl;
+                            if (this->partition.getBlock(predecessor) == block && isSilent(predecessor) && !stateLabels[this->partition.getPosition(predecessor) - block.getBeginIndex()].get(stateClassIndex)) {
+                                std::cout << "assigning label " << stateClassIndex << std::endl;
+                                stateStack.push_back(predecessor);
+                                stateLabels[this->partition.getPosition(predecessor) - block.getBeginIndex()].set(stateClassIndex);
+                            }
+                        }
+                    }
+                }
+            }
+            
+            return stateLabels;
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlockOfSplitterWeak(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            computeConditionalProbabilitiesForNonSilentStates(block);
+            // First, we need to compute a labeling of the states that expresses which of the non-silent blocks they can reach.
+            std::vector<uint_fast64_t> nonSilentBlockIndices = computeNonSilentBlocks(block);
+            std::cout << "indices: " << std::endl;
+            for (auto el : nonSilentBlockIndices) {
+                std::cout << el << std::endl;
+            }
+            
+            for (int ind = 0; ind < nonSilentBlockIndices.size() - 1; ++ind) {
+                for (int inner = nonSilentBlockIndices[ind]; inner < nonSilentBlockIndices[ind + 1]; ++inner) {
+                    std::cout << this->partition.getState(inner) << " is in class " << ind << std::endl;
+                }
+            }
+            
+            std::vector<storm::storage::BitVector> weakStateLabels = computeWeakStateLabelingBasedOnNonSilentBlocks(block, nonSilentBlockIndices);
+            
+            std::cout << "labels: " << std::endl;
+            for (auto el : weakStateLabels) {
+                std::cout << el << std::endl;
+            }
+            
+            // Then split the block according to this labeling.
+            // CAUTION: that this assumes that the positions of the states in the partition are not update until after
+            // the sorting is over. Otherwise, this interferes with the data used in the sorting process.
+            storm::storage::sparse::state_type originalBlockIndex = block.getBeginIndex();
+            auto result = this->partition.splitBlock(block,
+                                                     [&weakStateLabels,&block,originalBlockIndex,this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                                                         std::cout << "comparing states " << state1 << " and " << state2 << std::endl;
+                                                         std::cout << (this->partition.getPosition(state1) - originalBlockIndex) << " and " << (this->partition.getPosition(state2) - originalBlockIndex) << std::endl;
+                                                         std::cout << "size: " << weakStateLabels.size() << std::endl;
+                                                         std::cout << "begin is " << block.getBeginIndex() << std::endl;
+                                                         std::cout << "orig begin " << originalBlockIndex << std::endl;
+                                                         return weakStateLabels[this->partition.getPosition(state1) - originalBlockIndex] < weakStateLabels[this->partition.getPosition(state2) - originalBlockIndex];
+                                                     },
+                                                     [this, &splitterQueue] (bisimulation::Block<BlockDataType>& block) {
+                                                         updateSilentProbabilitiesBasedOnTransitions(block);
+                                                         
+                                                         // Insert the new block as a splitter.
+                                                         block.data().setSplitter();
+                                                         splitterQueue.emplace_back(&block);
+                                                     });
+            
+            // If the block was split, we also update the silent probabilities of the
+            if (result) {
+                updateSilentProbabilitiesBasedOnTransitions(block);
+                
+                // Insert the new block as a splitter.
+                block.data().setSplitter();
+                splitterQueue.emplace_back(&block);
+            } else {
+                block.resetMarkers();
+            }
+        }
+        
+        template<typename ModelType>
+        void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitterWeak(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            for (auto block : predecessorBlocks) {
+                std::cout << "found predecessor block " << block->getId() << std::endl;
+                if (*block != splitter) {
+                    refinePredecessorBlockOfSplitterWeak(*block, splitterQueue);
+                } else {
+                    // If the block to split is the splitter itself, we must not do any splitting here.
+                    block->resetMarkers();
+                }
+                
+                block->data().setNeedsRefinement(false);
+            }
+        }
+        
         template<typename ModelType>
         void DeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
             // The outline of the refinement is as follows.
@@ -297,10 +465,11 @@ namespace storm {
             // Finally, we use the information obtained in the first part for the actual splitting process in which all
             // predecessor blocks of the splitter are split based on the probabilities computed earlier.
             
-            // (1)
-            std::list<Block<BlockDataType>*> predecessorBlocks;
+            std::cout << "current partition is" << std::endl;
+            this->partition.print();
+            std::cout << "refining using splitter " << splitter.getId() << std::endl;
             
-            // (2)
+            std::list<Block<BlockDataType>*> predecessorBlocks;
             storm::storage::sparse::state_type currentPosition = splitter.getBeginIndex();
             bool splitterIsPredecessorBlock = false;
             for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.end(splitter); splitterIt != splitterIte; ++splitterIt, ++currentPosition) {
@@ -330,7 +499,7 @@ namespace storm {
                             splitterIsPredecessorBlock = true;
                             moveStateInSplitter(predecessor, predecessorBlock, currentPosition, elementsToSkip);
                         }
-
+                        
                         // Insert the block into the list of blocks to refine (if that has not already happened).
                         if (!predecessorBlock.data().needsRefinement()) {
                             predecessorBlocks.emplace_back(&predecessorBlock);
@@ -352,9 +521,15 @@ namespace storm {
             
             // Finally, we split the block based on the precomputed probabilities and the chosen bisimulation type.
             if (this->options.type == BisimulationType::Strong || this->model.getType() == storm::models::ModelType::Ctmc) {
-                refinePredecessorBlocksOfSplitter(predecessorBlocks, splitterQueue);
+                refinePredecessorBlocksOfSplitterStrong(predecessorBlocks, splitterQueue);
             } else {
-                assert(false);
+                // If the splitter is a predecessor of we can use the computed probabilities to update the silent
+                // probabilities.
+                if (splitterIsPredecessorBlock) {
+                    updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(splitter);
+                }
+                
+                refinePredecessorBlocksOfSplitterWeak(splitter, predecessorBlocks, splitterQueue);
             }
         }
         

src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h

@@ -42,7 +42,7 @@ namespace storm {
             virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) override;
 
         private:
-            virtual void refinePredecessorBlocksOfSplitter(std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            virtual void refinePredecessorBlocksOfSplitterStrong(std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
 
             /*!
              * Performs the necessary steps to compute a weak bisimulation on a DTMC.
@@ -69,6 +69,9 @@ namespace storm {
             // Retrieves whether the given state is silent.
             bool isSilent(storm::storage::sparse::state_type const& state) const;
             
+            // Retrieves whether the given state has a non-zero silent probability.
+            bool hasNonZeroSilentProbability(storm::storage::sparse::state_type const& state) const;
+            
             // Retrieves whether the given predecessor of the splitters possibly needs refinement.
             bool possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& predecessorBlock) const;
             
@@ -87,6 +90,29 @@ namespace storm {
             // Explores the remaining predecessors of the splitter.
             void exploreRemainingStatesOfSplitter(bisimulation::Block<BlockDataType>& splitter);
             
+            // Updates the silent probabilities of the states in the block based on the probabilities of going to the splitter.
+            void updateSilentProbabilitiesBasedOnProbabilitiesToSplitter(bisimulation::Block<BlockDataType>& block);
+
+            // Updates the silent probabilities of the states in the block based on a forward exploration of the transitions
+            // of the states.
+            void updateSilentProbabilitiesBasedOnTransitions(bisimulation::Block<BlockDataType>& block);
+            
+            // Refines the predecessor blocks of the splitter wrt. weak bisimulation in DTMCs.
+            void refinePredecessorBlocksOfSplitterWeak(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            // Refines the given block wrt to weak bisimulation in DTMCs.
+            void refinePredecessorBlockOfSplitterWeak(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            
+            // Converts the one-step probabilities of going into the splitter into the conditional probabilities needed
+            // for weak bisimulation (on DTMCs).
+            void computeConditionalProbabilitiesForNonSilentStates(bisimulation::Block<BlockDataType>& block);
+            
+            // Computes the (indices of the) blocks of non-silent states within the block.
+            std::vector<uint_fast64_t> computeNonSilentBlocks(bisimulation::Block<BlockDataType>& block);
+            
+            // Computes a labeling for all states of the block that identifies in which block they need to end up.
+            std::vector<storm::storage::BitVector> computeWeakStateLabelingBasedOnNonSilentBlocks(bisimulation::Block<BlockDataType> const& block, std::vector<uint_fast64_t> const& nonSilentBlockIndices);
+            
             // A vector that holds the probabilities of states going into the splitter. This is used by the method that
             // refines a block based on probabilities.
             std::vector<ValueType> probabilitiesToCurrentSplitter;

src/storage/bisimulation/Partition.cpp

@@ -107,13 +107,18 @@ namespace storm {
             }
             
             template<typename DataType>
-            void Partition<DataType>::mapStatesToPositions(Block<DataType> const& block) {
-                storm::storage::sparse::state_type position = block.getBeginIndex();
-                for (auto stateIt = this->begin(block), stateIte = this->end(block); stateIt != stateIte; ++stateIt, ++position) {
-                    this->positions[*stateIt] = position;
+            void Partition<DataType>::mapStatesToPositions(std::vector<storm::storage::sparse::state_type>::const_iterator first, std::vector<storm::storage::sparse::state_type>::const_iterator last) {
+                storm::storage::sparse::state_type position = std::distance(this->states.cbegin(), first);
+                for (; first != last; ++first, ++position) {
+                    this->positions[*first] = position;
                 }
             }
             
+            template<typename DataType>
+            void Partition<DataType>::mapStatesToPositions(Block<DataType> const& block) {
+                mapStatesToPositions(this->begin(block), this->end(block));
+            }
+            
             template<typename DataType>
             Block<DataType>& Partition<DataType>::getBlock(storm::storage::sparse::state_type state) {
                 return *this->stateToBlockMapping[state];
@@ -172,10 +177,42 @@ namespace storm {
                 return this->states.end();
             }
             
+            template<typename DataType>
+            void Partition<DataType>::sortRange(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions) {
+                std::sort(this->states.begin() + beginIndex, this->states.begin() + endIndex, less);
+                if (updatePositions) {
+                    mapStatesToPositions(this->states.begin() + beginIndex, this->states.begin() + endIndex);
+                }
+            }
+            
+            template<typename DataType>
+            void Partition<DataType>::sortBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions) {
+                sortRange(block.getBeginIndex(), block.getEndIndex(), less, updatePositions);
+            }
+            
+            template<typename DataType>
+            std::vector<uint_fast64_t> Partition<DataType>::computeRangesOfEqualValue(uint_fast64_t startIndex, uint_fast64_t endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less) {
+                auto it = this->states.cbegin() + startIndex;
+                auto ite = this->states.cbegin() + endIndex;
+                
+                std::vector<storm::storage::sparse::state_type>::const_iterator upperBound;
+                std::vector<uint_fast64_t> result;
+                result.push_back(startIndex);
+                do {
+                    upperBound = std::upper_bound(it, ite, *it, less);
+                    result.push_back(std::distance(this->states.cbegin(), upperBound));
+                    it = upperBound;
+                } while (upperBound != ite);
+                
+                return result;
+            }
+            
             template<typename DataType>
             std::pair<typename std::vector<std::unique_ptr<Block<DataType>>>::iterator, bool> Partition<DataType>::splitBlock(Block<DataType>& block, storm::storage::sparse::state_type position) {
                 STORM_LOG_THROW(position >= block.getBeginIndex() && position <= block.getEndIndex(), storm::exceptions::InvalidArgumentException, "Cannot split block at illegal position.");
 
+                std::cout << "splitting block " << block.getId() << " at pos " << position << std::endl;
+                
                 // In case one of the resulting blocks would be empty, we simply return the current block and do not create
                 // a new one.
                 if (position == block.getBeginIndex() || position == block.getEndIndex()) {
@@ -199,27 +236,38 @@ namespace storm {
             
             template<typename DataType>
             bool Partition<DataType>::splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback) {
-                // Sort the range of the block such that all states that have the label are moved to the front.
-                std::sort(this->begin(block), this->end(block), less);
+                // Sort the block, but leave the positions untouched.
+                this->sortBlock(block, less, false);
                 
-                // Update the positions vector.
-                mapStatesToPositions(block);
+                auto originalBegin = block.getBeginIndex();
+                auto originalEnd = block.getEndIndex();
+                
+                std::cout << "sorted block:" << std::endl;
+                for (auto stateIt = this->begin(block), stateIte = this->end(block); stateIt != stateIte; ++stateIt) {
+                    std::cout << *stateIt << std::endl;
+                }
+                
+                auto it = this->states.cbegin() + block.getBeginIndex();
+                auto ite = this->states.cbegin() + block.getEndIndex();
+                std::cout << "splitting between " << block.getBeginIndex() << " and " << block.getEndIndex() << std::endl;
                 
-                // Now we can check whether the block needs to be split, which is the case iff the changed function returns
-                // true for the first and last element of the remaining state range.
-                storm::storage::sparse::state_type begin = block.getBeginIndex();
-                storm::storage::sparse::state_type end = block.getEndIndex() - 1;
                 bool wasSplit = false;
-                while (less(states[begin], states[end])) {
-                    wasSplit = true;
-                    
-                    auto range = std::equal_range(states.begin() + begin, states.begin() + end, states[begin], less);
-                    begin = std::distance(states.begin(), range.second);
-                    auto result = this->splitBlock(block, begin);
-                    if (result.second) {
+                std::vector<storm::storage::sparse::state_type>::const_iterator upperBound;
+                do {
+                    std::cout << "it (" << *it << ") less than ite-1 (" << *(ite-1) << "? " << less(*it, *(ite - 1)) << std::endl;
+                    upperBound = std::upper_bound(it, ite, *it, less);
+                    std::cout << "upper bound is " << std::distance(this->states.cbegin(), upperBound);
+                    if (upperBound != ite) {
+                        wasSplit = true;
+                        auto result = this->splitBlock(block, std::distance(this->states.cbegin(), upperBound));
                         newBlockCallback(**result.first);
                     }
-                }
+                    it = upperBound;
+                } while (upperBound != ite);
+                
+                // Finally, repair the positions mapping.
+                mapStatesToPositions(this->states.begin() + originalBegin, this->states.begin() + originalEnd);
+
                 return wasSplit;
             }
             

src/storage/bisimulation/Partition.h

@@ -52,6 +52,15 @@ namespace storm {
                 // of the states.
                 std::pair<typename std::vector<std::unique_ptr<Block<DataType>>>::iterator, bool> splitBlock(Block<DataType>& block, storm::storage::sparse::state_type position);
 
+                // Sorts the given range of the partitition according to the given order.
+                void sortRange(storm::storage::sparse::state_type beginIndex, storm::storage::sparse::state_type endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions = true);
+                
+                // Sorts the block according to the given order.
+                void sortBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, bool updatePositions = true);
+
+                // Computes the start indices of equal ranges within the given range wrt. to the given less function.
+                std::vector<uint_fast64_t> computeRangesOfEqualValue(uint_fast64_t startIndex, uint_fast64_t endIndex, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less);
+                
                 // Splits the block by sorting the states according to the given function and then identifying the split
                 // points. The callback function is called for every newly created block.
                 bool splitBlock(Block<DataType>& block, std::function<bool (storm::storage::sparse::state_type, storm::storage::sparse::state_type)> const& less, std::function<void (Block<DataType>&)> const& newBlockCallback);
@@ -134,6 +143,9 @@ namespace storm {
                 
                 // Update the state to position for the states in the given block.
                 void mapStatesToPositions(Block<DataType> const& block);
+
+                // Update the state to position for the states in the given range.
+                void mapStatesToPositions(std::vector<storm::storage::sparse::state_type>::const_iterator first, std::vector<storm::storage::sparse::state_type>::const_iterator last);
                 
                 // Swaps the positions of the two states given by their positions.
                 void swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2);