tried different approach for bisim for MDPs

Former-commit-id: 92d56a4620

src/builder/DdPrismModelBuilder.cpp

@@ -1020,7 +1020,7 @@ namespace storm {
             
             preparedProgram = preparedProgram.substituteConstants();
             
-            STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl);
+            STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl);
             
             // Start by initializing the structure used for storing all information needed during the model generation.
             // In particular, this creates the meta variables used to encode the model.

src/builder/ExplicitPrismModelBuilder.cpp

@@ -295,7 +295,7 @@ namespace storm {
             // Now that the program is fixed, we we need to substitute all constants with their concrete value.
             preparedProgram = preparedProgram.substituteConstants();
             
-            STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl);
+            STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl);
                 
             // Select the appropriate reward models (after the constants have been substituted).
             std::vector<std::reference_wrapper<storm::prism::RewardModel const>> selectedRewardModels;

src/storage/Distribution.cpp

@@ -3,8 +3,14 @@
 #include <algorithm>
 #include <iostream>
 
+#include "src/utility/macros.h"
+#include "src/utility/constants.h"
+#include "src/utility/ConstantsComparator.h"
+
 #include "src/settings/SettingsManager.h"
 
+#include "src/adapters/CarlAdapter.h"
+
 namespace storm {
     namespace storage {
         
@@ -14,7 +20,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        bool Distribution<ValueType>::operator==(Distribution<ValueType> const& other) const {
+        bool Distribution<ValueType>::equals(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
             // We need to check equality by ourselves, because we need to account for epsilon differences.
             if (this->distribution.size() != other.distribution.size() || this->getHash() != other.getHash()) {
                 return false;
@@ -28,7 +34,7 @@ namespace storm {
                 if (first1->first != first2->first) {
                     return false;
                 }
-                if (std::abs(first1->second - first2->second) > 1e-6) {
+                if (!comparator.isEqual(first1->second, first2->second)) {
                     return false;
                 }
             }
@@ -38,10 +44,29 @@ namespace storm {
         
         template<typename ValueType>
         void Distribution<ValueType>::addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability) {
-            if (this->distribution.find(state) == this->distribution.end()) {
+            auto it = this->distribution.find(state);
+            if (it == this->distribution.end()) {
                 this->hash += static_cast<std::size_t>(state);
+                this->distribution.emplace_hint(it, state, probability);
+            } else {
+                it->second += probability;
+            }
+        }
+        
+        template<typename ValueType>
+        void Distribution<ValueType>::removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator) {
+            auto it = this->distribution.find(state);
+            STORM_LOG_ASSERT(it != this->distribution.end(), "Cannot remove probability, because the state is not in the support of the distribution.");
+            it->second -= probability;
+            if (comparator.isZero(it->second)) {
+                this->distribution.erase(it);
             }
-            this->distribution[state] += probability;
+        }
+        
+        template<typename ValueType>
+        void Distribution<ValueType>::shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator) {
+            removeProbability(fromState, probability, comparator);
+            addProbability(toState, probability);
         }
         
         template<typename ValueType>
@@ -68,7 +93,7 @@ namespace storm {
         void Distribution<ValueType>::scale(storm::storage::sparse::state_type const& state) {
             auto probabilityIterator = this->distribution.find(state);
             if (probabilityIterator != this->distribution.end()) {
-                ValueType scaleValue = 1 / probabilityIterator->second;
+                ValueType scaleValue = storm::utility::one<ValueType>() / probabilityIterator->second;
                 this->distribution.erase(probabilityIterator);
                 
                 for (auto& entry : this->distribution) {
@@ -82,6 +107,11 @@ namespace storm {
             return this->hash ^ (this->distribution.size() << 8);
         }
         
+        template<typename ValueType>
+        std::size_t Distribution<ValueType>::size() const {
+            return this->distribution.size();
+        }
+        
         template<typename ValueType>
         std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution) {
             out << "{";
@@ -93,7 +123,34 @@ namespace storm {
             return out;
         }
         
+        template<typename ValueType>
+        bool Distribution<ValueType>::less(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
+            if (this->size() != other.size()) {
+                return this->size() < other.size();
+            }
+            auto firstIt = this->begin();
+            auto firstIte = this->end();
+            auto secondIt = other.begin();
+            for (; firstIt != firstIte; ++firstIt, ++secondIt) {
+                // If the two blocks already differ, we can decide which distribution is smaller.
+                if (firstIt->first != secondIt->first) {
+                    return firstIt->first < secondIt->first;
+                }
+                
+                // If the blocks are the same, but the probability differs, we can also decide which distribution is smaller.
+                if (!comparator.isEqual(firstIt->second, secondIt->second)) {
+                    return comparator.isLess(firstIt->second, secondIt->second);
+                }
+            }
+            return false;
+        }
+        
         template class Distribution<double>;
         template std::ostream& operator<<(std::ostream& out, Distribution<double> const& distribution);
+        
+#ifdef STORM_HAVE_CARL
+        template class Distribution<storm::RationalFunction>;
+        template std::ostream& operator<<(std::ostream& out, Distribution<storm::RationalFunction> const& distribution);
+#endif
     }
 }

src/storage/Distribution.h

@@ -8,6 +8,11 @@
 #include "src/storage/sparse/StateType.h"
 
 namespace storm {
+    namespace utility {
+        template <typename ValueType>
+        class ConstantsComparator;
+    }
+    
     namespace storage {
         
         template<typename ValueType>
@@ -28,7 +33,7 @@ namespace storm {
              * @param other The distribution with which the current distribution is to be compared.
              * @return True iff the two distributions are equal.
              */
-            bool operator==(Distribution const& other) const;
+            bool equals(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>()) const;
             
             /*!
              * Assigns the given state the given probability under this distribution.
@@ -38,13 +43,34 @@ namespace storm {
              */
             void addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability);
             
+            /*!
+             * Removes the given probability mass of going to the given state.
+             *
+             * @param state The state for which to remove the probability.
+             * @param probability The probability to remove.
+             * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the
+             * entry is removed.
+             */
+            void removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
+            
+            /*!
+             * Removes the probability mass from one state and adds it to another.
+             *
+             * @param fromState The state from which to take the probability mass.
+             * @param toState The state from which to which to add the probability mass.
+             * @param probability The probability mass to shift.
+             * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the
+             * entry is removed.
+             */
+            void shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
+            
             /*!
              * Retrieves an iterator to the elements in this distribution.
              *
              * @return The iterator to the elements in this distribution.
              */
             iterator begin();
-
+            
             /*!
              * Retrieves an iterator to the elements in this distribution.
              *
@@ -58,7 +84,7 @@ namespace storm {
              * @return The iterator past the elements in this distribution.
              */
             iterator end();
-
+            
             /*!
              * Retrieves an iterator past the elements in this distribution.
              *
@@ -82,16 +108,23 @@ namespace storm {
              */
             std::size_t getHash() const;
             
+            /*!
+             * Retrieves the size of the distribution, i.e. the size of the support set.
+             */
+            std::size_t size() const;
+            
+            bool less(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const;
+            
         private:
             // A list of states and the probabilities that are assigned to them.
             container_type distribution;
             
-            // A hash value that is maintained to allow for quicker equality comparison between distribution.s
+            // A hash value that is maintained to allow for quicker equality comparison between distributions.
             std::size_t hash;
         };
         
         template<typename ValueType>
-            std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution);
+        std::ostream& operator<<(std::ostream& out, Distribution<ValueType> const& distribution);
     }
 }
 

src/storage/bisimulation/BisimulationDecomposition.cpp

@@ -4,6 +4,7 @@
 
 #include "src/models/sparse/Dtmc.h"
 #include "src/models/sparse/Ctmc.h"
+#include "src/models/sparse/Mdp.h"
 #include "src/models/sparse/StandardRewardModel.h"
 
 #include "src/storage/bisimulation/DeterministicBlockData.h"
@@ -167,8 +168,12 @@ namespace storm {
             } else {
                 this->initializeLabelBasedPartition();
             }
+            std::cout << "initial partition is " << std::endl;
+            this->partition.print();
             std::chrono::high_resolution_clock::duration initialPartitionTime = std::chrono::high_resolution_clock::now() - initialPartitionStart;
             
+            this->initialize();
+            
             std::chrono::high_resolution_clock::time_point refinementStart = std::chrono::high_resolution_clock::now();
             this->performPartitionRefinement();
             std::chrono::high_resolution_clock::duration refinementTime = std::chrono::high_resolution_clock::now() - refinementStart;
@@ -274,6 +279,11 @@ namespace storm {
             }
         }
         
+        template<typename ModelType, typename BlockDataType>
+        void BisimulationDecomposition<ModelType, BlockDataType>::initialize() {
+            // Intentionally left empty.
+        }
+        
         template<typename ModelType, typename BlockDataType>
         void BisimulationDecomposition<ModelType, BlockDataType>::extractDecompositionBlocks() {
             // Now move the states from the internal partition into their final place in the decomposition. We do so in
@@ -290,10 +300,12 @@ namespace storm {
         
         template class BisimulationDecomposition<storm::models::sparse::Dtmc<double>, bisimulation::DeterministicBlockData>;
         template class BisimulationDecomposition<storm::models::sparse::Ctmc<double>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Mdp<double>, bisimulation::DeterministicBlockData>;
         
 #ifdef STORM_HAVE_CARL
         template class BisimulationDecomposition<storm::models::sparse::Dtmc<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
         template class BisimulationDecomposition<storm::models::sparse::Ctmc<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
+        template class BisimulationDecomposition<storm::models::sparse::Mdp<storm::RationalFunction>, bisimulation::DeterministicBlockData>;
 #endif
     }
 }

src/storage/bisimulation/BisimulationDecomposition.h

@@ -185,6 +185,11 @@ namespace storm {
              */
             virtual void initializeMeasureDrivenPartition();
             
+            /*!
+             * A function that can initialize auxiliary data structures. It is called after initializing the initial partition.
+             */
+            virtual void initialize();
+            
             /*!
              * Computes the set of states with probability 0/1 for satisfying phi until psi. This is used for the measure
              * driven initial partition.

src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp

@@ -172,7 +172,7 @@ namespace storm {
                 
                 split |= this->partition.splitBlock(*blockToRefineProbabilistically,
                                                     [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
-                                                        return getProbabilityToSplitter(state1) < getProbabilityToSplitter(state2);
+                                                        return this->comparator.isLess(getProbabilityToSplitter(state1), getProbabilityToSplitter(state2));
                                                     },
                                                     [&splitterQueue] (Block<BlockDataType>& block) {
                                                         splitterQueue.emplace_back(&block); block.data().setSplitter();

src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp

@@ -5,19 +5,21 @@
 
 #include "src/utility/graph.h"
 
+#include "src/utility/macros.h"
 #include "src/exceptions/IllegalFunctionCallException.h"
 
+#include "src/adapters/CarlAdapter.h"
+
 namespace storm {
     namespace storage {
-
+        
         using namespace bisimulation;
         
         template<typename ModelType>
-        NondeterministicModelBisimulationDecomposition<ModelType>::NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>::Options const& options) : BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>(model, model.getTransitionMatrix().transpose(false), options), choiceToStateMapping(model.getNumberOfChoices()), probabilitiesToCurrentSplitter(model.getNumberOfChoices(), storm::utility::zero<ValueType>()) {
+        NondeterministicModelBisimulationDecomposition<ModelType>::NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>::Options const& options) : BisimulationDecomposition<ModelType, NondeterministicModelBisimulationDecomposition::BlockDataType>(model, model.getTransitionMatrix().transpose(false), options), choiceToStateMapping(model.getNumberOfChoices()), quotientDistributions(model.getNumberOfChoices()), orderedQuotientDistributions(model.getNumberOfChoices()) {
             STORM_LOG_THROW(options.type == BisimulationType::Strong, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation is currently not supported for nondeterministic models.");
-            this->createChoiceToStateMapping();
         }
-    
+        
         template<typename ModelType>
         std::pair<storm::storage::BitVector, storm::storage::BitVector> NondeterministicModelBisimulationDecomposition<ModelType>::getStatesWithProbability01() {
             STORM_LOG_THROW(static_cast<bool>(this->options.optimalityType), storm::exceptions::IllegalFunctionCallException, "Can only compute states with probability 0/1 with an optimization direction (min/max).");
@@ -27,7 +29,13 @@ namespace storm {
                 return storm::utility::graph::performProb01Max(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(), this->model.getBackwardTransitions(), this->options.phiStates.get(), this->options.psiStates.get());
             }
         }
-
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::initialize() {
+            this->createChoiceToStateMapping();
+            this->initializeQuotientDistributions();
+        }
+        
         template<typename ModelType>
         void NondeterministicModelBisimulationDecomposition<ModelType>::createChoiceToStateMapping() {
             std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
@@ -39,192 +47,172 @@ namespace storm {
         }
         
         template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
-            STORM_LOG_ASSERT(false, "Not yet implemented");
-        }
-        
-        template<typename ModelType>
-        bool NondeterministicModelBisimulationDecomposition<ModelType>::possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& predecessorBlock) const {
-            return predecessorBlock.getNumberOfStates() > 1 && !predecessorBlock.data().absorbing();
-        }
-
-        template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::clearProbabilitiesToSplitter(storm::storage::sparse::state_type state) {
+        void NondeterministicModelBisimulationDecomposition<ModelType>::initializeQuotientDistributions() {
             std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
-            for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
-                probabilitiesToCurrentSplitter[choice] = storm::utility::zero<ValueType>();
+            for (auto choice = 0; choice < nondeterministicChoiceIndices.back(); ++choice) {
+                for (auto entry : this->model.getTransitionMatrix().getRow(choice)) {
+                    if (!this->comparator.isZero(entry.getValue())) {
+                        this->quotientDistributions[choice].addProbability(this->partition.getBlock(entry.getColumn()).getId(), entry.getValue());
+                    }
+                }
+                orderedQuotientDistributions[choice] = &this->quotientDistributions[choice];
             }
-        }
-        
-        template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::increaseProbabilityToSplitter(storm::storage::sparse::state_type state, uint_fast64_t choice, bisimulation::Block<BlockDataType> const& predecessorBlock, ValueType const& value) {
-            storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(state);
             
-            // If the position of the state is to the right of marker1, we have not seen it before. That means we need to
-            // clear the probability associated with all choices of the given state.
-            if (predecessorPosition >= predecessorBlock.data().marker1()) {
-                clearProbabilitiesToSplitter(state);
+            for (auto state = 0; state < this->model.getNumberOfStates(); ++state) {
+                updateOrderedQuotientDistributions(state);
             }
-            
-            // Now increase the probability of the given choice of the given state.
-            probabilitiesToCurrentSplitter[choice] += value;
-        }
-        
-        template <typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock) {
-            this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1()));
-            predecessorBlock.data().incrementMarker1();
         }
         
-        template <typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock) {
-            this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker2()));
-            predecessorBlock.data().incrementMarker2();
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::updateOrderedQuotientDistributions(storm::storage::sparse::state_type state) {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            std::sort(this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state], this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state + 1],
+                      [this] (storm::storage::Distribution<ValueType> const* dist1, storm::storage::Distribution<ValueType> const* dist2) {
+                          return dist1->less(*dist2, this->comparator);
+                      });
         }
         
-        template <typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip) {
-            storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
-            
-            // If the predecessors of the given predecessor were already explored, we can move it easily.
-            if (predecessorPosition <= currentPositionInSplitter + elementsToSkip) {
-                this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1()));
-                predecessorBlock.data().incrementMarker1();
-            } else {
-                // Otherwise, we need to move the predecessor, but we need to make sure that we explore its
-                // predecessors later. We do this by moving it to a range at the beginning of the block that will hold
-                // all predecessors in the splitter whose predecessors have yet to be explored.
-                if (predecessorBlock.data().marker2() == predecessorBlock.data().marker1()) {
-                    this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition);
-                    this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1);
-                } else {
-                    this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition);
-                    this->partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.data().marker1());
-                    this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1);
-                }
-                
-                // Since we had to move an already explored state to the right of the current position,
-                ++elementsToSkip;
-                predecessorBlock.data().incrementMarker1();
-                predecessorBlock.data().incrementMarker2();
-            }
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
+            std::cout << "Found " << this->partition.size() << " blocks" << std::endl;
+            this->partition.print();
+            STORM_LOG_ASSERT(false, "Not yet implemented");
         }
         
         template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::insertIntoPredecessorList(bisimulation::Block<BlockDataType>& predecessorBlock, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks) {
-            // Insert the block into the list of blocks to refine (if that has not already happened).
-            if (!predecessorBlock.data().needsRefinement()) {
-                predecessorBlocks.emplace_back(&predecessorBlock);
-                predecessorBlock.data().setNeedsRefinement();
-            }
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& block) const {
+            return block.getNumberOfStates() > 1 && !block.data().absorbing();
         }
         
-        template <typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::exploreRemainingStatesOfSplitter(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks) {
-            for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.begin(splitter) + (splitter.data().marker2() - splitter.getBeginIndex()); splitterIt != splitterIte; ++splitterIt) {
-                storm::storage::sparse::state_type currentState = *splitterIt;
-                
-                for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
-                    storm::storage::sparse::state_type choice = predecessorEntry.getColumn();
-                    storm::storage::sparse::state_type predecessor = choiceToStateMapping[choice];
-                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessor);
-                    
-                    // If the block does not need to be refined, we skip it.
-                    if (!possiblyNeedsRefinement(predecessorBlock)) {
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::updateQuotientDistributionsOfPredecessors(Block<BlockDataType> const& newBlock, Block<BlockDataType> const& oldBlock, std::deque<Block<BlockDataType>*>& splitterQueue) {
+            uint_fast64_t lastState = 0;
+            bool lastStateInitialized = false;
+            
+            for (auto stateIt = this->partition.begin(newBlock), stateIte = this->partition.end(newBlock); stateIt != stateIte; ++stateIt) {
+                for (auto predecessorEntry : this->backwardTransitions.getRow(*stateIt)) {
+                    if (this->comparator.isZero(predecessorEntry.getValue())) {
                         continue;
                     }
                     
-                    // We keep track of the probability of the predecessor moving to the splitter.
-                    increaseProbabilityToSplitter(predecessor, choice, predecessorBlock, predecessorEntry.getValue());
+                    storm::storage::sparse::state_type predecessorChoice = predecessorEntry.getColumn();
+                    storm::storage::sparse::state_type predecessorState = choiceToStateMapping[predecessorChoice];
+                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessorState);
+                    
+                    // If the predecessor block is not marked as to-refined, we do so now.
+                    if (!predecessorBlock.data().splitter()) {
+                        predecessorBlock.data().setSplitter();
+                        splitterQueue.push_back(&predecessorBlock);
+                    }
                     
-                    // Only move the state if it has not been seen as a predecessor before.
-                    storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
-                    if (predecessorPosition >= predecessorBlock.data().marker1()) {
-                        moveStateToMarker1(predecessor, predecessorBlock);
+                    if (lastStateInitialized) {
+                        // If we have skipped to the choices of the next state, we need to repair the order of the
+                        // distributions for the last state.
+                        if (lastState != predecessorState) {
+                            updateOrderedQuotientDistributions(lastState);
+                            lastState = predecessorState;
+                        }
+                    } else {
+                        lastStateInitialized = true;
+                        lastState = choiceToStateMapping[predecessorChoice];
                     }
                     
-                    insertIntoPredecessorList(predecessorBlock, predecessorBlocks);
+                    // Now shift the probability from this transition from the old block to the new one.
+                    this->quotientDistributions[predecessorChoice].shiftProbability(oldBlock.getId(), newBlock.getId(), predecessorEntry.getValue());
                 }
             }
             
-            // Finally, we can reset the second marker.
-            splitter.data().setMarker2(splitter.getBeginIndex());
+            if (lastStateInitialized) {
+                updateOrderedQuotientDistributions(lastState);
+            }
         }
         
         template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitter(std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
-            // TODO
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkQuotientDistributions() const {
+            
         }
         
         template<typename ModelType>
-        void NondeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
-            // The outline of the refinement is as follows.
-            //
-            // We iterate over all states of the splitter and determine for each predecessor the state the probability
-            // entering the splitter. These probabilities are written to a member vector so that after the iteration
-            // process we have the probabilities of all predecessors of the splitter of entering the splitter in one
-            // step. To directly separate the states having a transition into the splitter from the ones who do not,
-            // we move the states to certain locations. That is, on encountering a predecessor of the splitter, it is
-            // moved to the beginning of its block. If the predecessor is in the splitter itself, we have to be a bit
-            // careful about where to move states.
-            //
-            // After this iteration, there may be states of the splitter whose predecessors have not yet been explored,
-            // so this needs to be done now.
-            //
-            // Finally, we use the information obtained in the first part for the actual splitting process in which all
-            // predecessor blocks of the splitter are split based on the probabilities computed earlier.
-            std::list<Block<BlockDataType>*> predecessorBlocks;
-            storm::storage::sparse::state_type currentPosition = splitter.getBeginIndex();
-            bool splitterIsPredecessorBlock = false;
-            for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.end(splitter); splitterIt != splitterIte; ++splitterIt, ++currentPosition) {
-                storm::storage::sparse::state_type currentState = *splitterIt;
-                
-                uint_fast64_t elementsToSkip = 0;
-                for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) {
-                    storm::storage::sparse::state_type choice = predecessorEntry.getColumn();
-                    storm::storage::sparse::state_type predecessor = choiceToStateMapping[choice];
-                    storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor);
-                    Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessor);
-                    
-                    // If the block does not need to be refined, we skip it.
-                    if (!possiblyNeedsRefinement(predecessorBlock)) {
-                        continue;
-                    }
-                    
-                    // We keep track of the probability of the predecessor moving to the splitter.
-                    increaseProbabilityToSplitter(predecessor, choice, predecessorBlock, predecessorEntry.getValue());
-                    
-                    // We only need to move the predecessor if its not already known as a predecessor already.
-                    if (predecessorPosition >= predecessorBlock.data().marker1()) {
-                        // If the predecessor block is not the splitter, we can move the state easily.
-                        if (predecessorBlock != splitter) {
-                            moveStateToMarker1(predecessor, predecessorBlock);
-                        } else {
-                            // If the predecessor is in the splitter, we need to be a bit more careful.
-                            splitterIsPredecessorBlock = true;
-                            moveStateInSplitter(predecessor, predecessorBlock, currentPosition, elementsToSkip);
-                        }
-                        
-                        insertIntoPredecessorList(predecessorBlock, predecessorBlocks);
-                    }
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::splitBlockAccordingToCurrentQuotientDistributions(Block<BlockDataType>& block, std::deque<Block<BlockDataType>*>& splitterQueue) {
+            bool split = this->partition.splitBlock(block,
+                                                    [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+                                                        return quotientDistributionsLess(state1, state2);
+                                                    },
+                                                    [this, &block, &splitterQueue] (Block<BlockDataType>& newBlock) {
+                                                        updateQuotientDistributionsOfPredecessors(newBlock, block, splitterQueue);
+                                                    });
+            
+            // The quotient distributions of the predecessors of block do not need to be updated, since the probability
+            // will go to the block with the same id as before.
+            
+//            std::cout << "partition after split: " << std::endl;
+//            this->partition.print();
+            
+            this->checkQuotientDistributions();
+            
+            return split;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+            STORM_LOG_TRACE("Comparing the quotient distributions of state " << state1 << " and " << state2 << ".");
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            
+//            std::cout << "state " << state1 << std::endl;
+//            for (int c = nondeterministicChoiceIndices[state1]; c < nondeterministicChoiceIndices[state1 + 1]; ++c) {
+//                std::cout << quotientDistributions[c] << std::endl;
+//            }
+//
+//            std::cout << "state " << state2 << std::endl;
+//            for (int c = nondeterministicChoiceIndices[state2]; c < nondeterministicChoiceIndices[state2 + 1]; ++c) {
+//                std::cout << quotientDistributions[c] << std::endl;
+//            }
+            
+            auto firstIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1];
+            auto firstIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1 + 1];
+            auto secondIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2];
+            auto secondIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2 + 1];
+            
+            for (; firstIt != firstIte && secondIt != secondIte; ++firstIt, ++secondIt) {
+                // If the current distributions are in a less-than relationship, we can return a result.
+                if ((*firstIt)->less(**secondIt, this->comparator)) {
+                    return true;
+                } else if ((*secondIt)->less(**firstIt, this->comparator)) {
+                    return false;
                 }
                 
-                // If, as a consequence of shifting states, we need to skip some elements, do so now.
-                splitterIt += elementsToSkip;
-                currentPosition += elementsToSkip;
+                // If the distributions matched, we need to advance both distribution iterators to the next distribution
+                // that is larger.
+                while (firstIt != firstIte && std::next(firstIt) != firstIte && !(*firstIt)->less(**std::next(firstIt), this->comparator)) {
+                    ++firstIt;
+                }
+                while (secondIt != secondIte && std::next(secondIt) != secondIte && !(*secondIt)->less(**std::next(secondIt), this->comparator)) {
+                    ++secondIt;
+                }
             }
             
-            // If the splitter was a predecessor block of itself, we potentially need to explore some states that have
-            // not been explored yet.
-            if (splitterIsPredecessorBlock) {
-                exploreRemainingStatesOfSplitter(splitter, predecessorBlocks);
+            if (firstIt == firstIte && secondIt != secondIte) {
+                return true;
+            }
+            return false;
+        }
+        
+        template<typename ModelType>
+        void NondeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            if (!possiblyNeedsRefinement(splitter)) {
+                return;
             }
             
-            // Finally, we split the block based on the precomputed probabilities and the chosen bisimulation type.
-            refinePredecessorBlocksOfSplitter(predecessorBlocks, splitterQueue);
+            STORM_LOG_TRACE("Refining block " << splitter.getId());
+            
+            splitBlockAccordingToCurrentQuotientDistributions(splitter, splitterQueue);
         }
         
         template class NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>;
         
+#ifdef STORM_HAVE_CARL
+        template class NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<storm::RationalFunction>>;
+#endif
+        
     }
 }

src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h

@@ -4,6 +4,8 @@
 #include "src/storage/bisimulation/BisimulationDecomposition.h"
 #include "src/storage/bisimulation/DeterministicBlockData.h"
 
+#include "src/storage/Distribution.h"
+
 namespace storm {
     namespace utility {
         template <typename ValueType> class ConstantsComparator;
@@ -37,44 +39,41 @@ namespace storm {
             
             virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) override;
             
+            virtual void initialize() override;
+            
         private:
             // Creates the mapping from the choice indices to the states.
             void createChoiceToStateMapping();
             
-            // Retrieves whether the given predecessor of the splitters possibly needs refinement.
-            bool possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& predecessorBlock) const;
-            
-            // Increases the probability of moving to the current splitter for the given choice of the given state.
-            void increaseProbabilityToSplitter(storm::storage::sparse::state_type state, uint_fast64_t choice, bisimulation::Block<BlockDataType> const& predecessorBlock, ValueType const& value);
-            
-            // Clears the probabilities of all choices of the given state.
-            void clearProbabilitiesToSplitter(storm::storage::sparse::state_type state);
+            // Initializes the quotient distributions wrt. to the current partition.
+            void initializeQuotientDistributions();
             
-            // Moves the given state to the position marked by marker1 moves the marker one step further.
-            void moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock);
+            // Retrieves whether the given block possibly needs refinement.
+            bool possiblyNeedsRefinement(bisimulation::Block<BlockDataType> const& block) const;
             
-            // Moves the given state to the position marked by marker2 the marker one step further.
-            void moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock);
+            // Splits the given block according to the current quotient distributions.
+            bool splitBlockAccordingToCurrentQuotientDistributions(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
             
-            // Moves the given state to a proper place in the splitter, depending on where the predecessor is located.
-            void moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block<BlockDataType>& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip);
+            // Retrieves whether the quotient distributions of state 1 are considered to be less than the ones of state 2.
+            bool quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2);
             
-            // Inserts the block into the list of predecessors if it is not already contained.
-            void insertIntoPredecessorList(bisimulation::Block<BlockDataType>& predecessorBlock, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks);
+            // Updates the ordered list of quotient distribution for the given state.
+            void updateOrderedQuotientDistributions(storm::storage::sparse::state_type state);
             
-            // Explores the remaining states of the splitter.
-            void exploreRemainingStatesOfSplitter(bisimulation::Block<BlockDataType>& splitter, std::list<bisimulation::Block<BlockDataType>*>& predecessorBlocks);
+            // Updates the quotient distributions of the predecessors of the new block by taking the probability mass
+            // away from the old block.
+            void updateQuotientDistributionsOfPredecessors(bisimulation::Block<BlockDataType> const& newBlock, bisimulation::Block<BlockDataType> const& oldBlock, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
             
-            // Refines the predecessor blocks wrt. strong bisimulation.
-            void refinePredecessorBlocksOfSplitter(std::list<bisimulation::Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
+            bool checkQuotientDistributions() const;
             
             // A mapping from choice indices to the state state that has this choice.
             std::vector<storm::storage::sparse::state_type> choiceToStateMapping;
             
-            // A vector that holds the probabilities for all nondeterministic choices of all states of going into the
-            // splitter. This is used by the method that refines a block based on probabilities.
-            std::vector<ValueType> probabilitiesToCurrentSplitter;
+            // A vector that holds the quotient distributions for all nondeterministic choices of all states.
+            std::vector<storm::storage::Distribution<ValueType>> quotientDistributions;
             
+            // A vector that stores for each state the ordered list of quotient distributions.
+            std::vector<storm::storage::Distribution<ValueType> const*> orderedQuotientDistributions;
         };
     }
 }

src/utility/ConstantsComparator.cpp

@@ -35,6 +35,11 @@ namespace storm {
             return false;
         }
         
+        template<typename ValueType>
+        bool ConstantsComparator<ValueType>::isLess(ValueType const& value1, ValueType const& value2) const {
+            return value1 < value2;
+        }
+        
         ConstantsComparator<float>::ConstantsComparator() : precision(static_cast<float>(storm::settings::generalSettings().getPrecision())) {
             // Intentionally left empty.
         }
@@ -63,6 +68,10 @@ namespace storm {
             return value == storm::utility::infinity<float>();
         }
         
+        bool ConstantsComparator<float>::isLess(float const& value1, float const& value2) const {
+            return std::abs(value1 - value2) < precision;
+        }
+        
         ConstantsComparator<double>::ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) {
             // Intentionally left empty.
         }
@@ -90,7 +99,11 @@ namespace storm {
         bool ConstantsComparator<double>::isConstant(double const& value) const {
             return true;
         }
-                
+        
+        bool ConstantsComparator<double>::isLess(double const& value1, double const& value2) const {
+            return std::abs(value1 - value2) < precision;
+        }
+        
         // Explicit instantiations.
         template class ConstantsComparator<double>;
         template class ConstantsComparator<float>;

src/utility/ConstantsComparator.h

@@ -21,6 +21,8 @@ namespace storm {
             bool isConstant(ValueType const& value) const;
             
             bool isInfinity(ValueType const& value) const;
+            
+            bool isLess(ValueType const& value1, ValueType const& value2) const;
         };
         
         // For floats we specialize this class and consider the comparison modulo some predefined precision.
@@ -41,6 +43,8 @@ namespace storm {
             
             bool isInfinity(float const& value) const;
 
+            bool isLess(float const& value1, float const& value2) const;
+            
         private:
             // The precision used for comparisons.
             float precision;
@@ -64,6 +68,8 @@ namespace storm {
             
             bool isConstant(double const& value) const;
             
+            bool isLess(double const& value1, double const& value2) const;
+            
         private:
             // The precision used for comparisons.
             double precision;

src/utility/storm.h

@@ -45,6 +45,7 @@
 
 // Headers for model processing.
 #include "src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h"
+#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
 
 // Headers for model checking.
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
@@ -114,7 +115,7 @@ namespace storm {
     }
     
     template<typename ModelType>
-    std::shared_ptr<ModelType> performSparseBisimulationMinimization(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+    std::shared_ptr<ModelType> performDeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
         std::cout << "Performing bisimulation minimization... ";
         typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options options;
         if (!formulas.empty()) {
@@ -132,19 +133,40 @@ namespace storm {
         return model;
     }
     
+    template<typename ModelType>
+    std::shared_ptr<ModelType> performNondeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+        std::cout << "Performing bisimulation minimization... ";
+        typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options options;
+        if (!formulas.empty()) {
+            options = typename storm::storage::NondeterministicModelBisimulationDecomposition<ModelType>::Options(*model, formulas);
+        }
+        if (storm::settings::bisimulationSettings().isWeakBisimulationSet()) {
+            options.type = storm::storage::BisimulationType::Weak;
+            options.bounded = false;
+        }
+        
+        storm::storage::NondeterministicModelBisimulationDecomposition<ModelType> bisimulationDecomposition(*model, options);
+        bisimulationDecomposition.computeBisimulationDecomposition();
+        model = bisimulationDecomposition.getQuotient();
+        std::cout << "done." << std::endl << std::endl;
+        return model;
+    }
+    
     template<typename ModelType, typename ValueType = typename ModelType::ValueType, typename std::enable_if<std::is_base_of<storm::models::sparse::Model<ValueType>, ModelType>::value, bool>::type = 0>
     std::shared_ptr<storm::models::ModelBase> preprocessModel(std::shared_ptr<ModelType> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
         if (storm::settings::generalSettings().isBisimulationSet()) {
             
             STORM_LOG_THROW(model->isSparseModel(), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for sparse models.");
-            STORM_LOG_THROW(model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs and CTMCs.");
+            STORM_LOG_THROW(model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc) || model->isOfType(storm::models::ModelType::Mdp), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs, CTMCs and MDPs.");
 
             model->reduceToStateBasedRewards();
 
             if (model->isOfType(storm::models::ModelType::Dtmc)) {
-                return performSparseBisimulationMinimization<storm::models::sparse::Dtmc<ValueType>>(model->template as<storm::models::sparse::Dtmc<ValueType>>(), formulas);
+                return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Dtmc<ValueType>>(model->template as<storm::models::sparse::Dtmc<ValueType>>(), formulas);
+            } else if (model->isOfType(storm::models::ModelType::Ctmc)) {
+                return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas);
             } else {
-                return performSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas);
+                return performNondeterministicSparseBisimulationMinimization<storm::models::sparse::Mdp<ValueType>>(model->template as<storm::models::sparse::Mdp<ValueType>>(), formulas);
             }
             
         }