Explicit MDP bisim working but unfortunately slow :(

Former-commit-id: 6714bdbd61

src/storage/Distribution.cpp

@@ -15,14 +15,14 @@ namespace storm {
     namespace storage {
         
         template<typename ValueType>
-        Distribution<ValueType>::Distribution() : hash(0) {
+        Distribution<ValueType>::Distribution() {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
         bool Distribution<ValueType>::equals(Distribution<ValueType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const {
             // We need to check equality by ourselves, because we need to account for epsilon differences.
-            if (this->distribution.size() != other.distribution.size() || this->getHash() != other.getHash()) {
+            if (this->distribution.size() != other.distribution.size()) {
                 return false;
             }
             
@@ -32,9 +32,11 @@ namespace storm {
             
             for (; first1 != last1; ++first1, ++first2) {
                 if (first1->first != first2->first) {
+                    std::cout << "false in first" << std::endl;
                     return false;
                 }
                 if (!comparator.isEqual(first1->second, first2->second)) {
+                    std::cout << "false in second " << std::setprecision(15) << first1->second << " vs " << first2->second << std::endl;
                     return false;
                 }
             }
@@ -46,7 +48,6 @@ namespace storm {
         void Distribution<ValueType>::addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability) {
             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;
@@ -101,12 +102,7 @@ namespace storm {
                 }
             }
         }
-        
-        template<typename ValueType>
-        std::size_t Distribution<ValueType>::getHash() const {
-            return this->hash ^ (this->distribution.size() << 8);
-        }
-        
+                
         template<typename ValueType>
         std::size_t Distribution<ValueType>::size() const {
             return this->distribution.size();

src/storage/Distribution.h

@@ -5,6 +5,8 @@
 #include <ostream>
 #include <boost/container/flat_map.hpp>
 
+#include "src/utility/OsDetection.h"
+
 #include "src/storage/sparse/StateType.h"
 
 namespace storm {
@@ -26,6 +28,14 @@ namespace storm {
              * Creates an empty distribution.
              */
             Distribution();
+
+            Distribution(Distribution const& other) = default;
+            Distribution& operator=(Distribution const& other) = default;
+
+#ifndef WINDOWS
+            Distribution(Distribution&& other) = default;
+            Distribution& operator=(Distribution&& other) = default;
+#endif
             
             /*!
              * Checks whether the two distributions specify the same probabilities to go to the same states.
@@ -101,13 +111,6 @@ namespace storm {
              */
             void scale(storm::storage::sparse::state_type const& state);
             
-            /*!
-             * Retrieves the hash value of the distribution.
-             *
-             * @return The hash value of the distribution.
-             */
-            std::size_t getHash() const;
-            
             /*!
              * Retrieves the size of the distribution, i.e. the size of the support set.
              */
@@ -118,9 +121,6 @@ namespace storm {
         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 distributions.
-            std::size_t hash;
         };
         
         template<typename ValueType>

src/storage/bisimulation/BisimulationDecomposition.cpp

@@ -84,7 +84,9 @@ namespace storm {
             std::shared_ptr<storm::logic::Formula const> newFormula = formula.asSharedPointer();
             
             if (formula.isProbabilityOperatorFormula()) {
-                optimalityType = formula.asProbabilityOperatorFormula().getOptimalityType();
+                if (formula.asProbabilityOperatorFormula().hasOptimalityType()) {
+                    optimalityType = formula.asProbabilityOperatorFormula().getOptimalityType();
+                }
                 newFormula = formula.asProbabilityOperatorFormula().getSubformula().asSharedPointer();
             } else if (formula.isRewardOperatorFormula()) {
                 optimalityType = formula.asRewardOperatorFormula().getOptimalityType();
@@ -168,8 +170,6 @@ 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();
@@ -218,9 +218,7 @@ namespace storm {
             uint_fast64_t iterations = 0;
             while (!splitterQueue.empty()) {
                 ++iterations;
-                // Optionally: sort the splitter queue according to some criterion (here: prefer small splitters).
-                std::sort(splitterQueue.begin(), splitterQueue.end(), [] (Block<BlockDataType> const* b1, Block<BlockDataType> const* b2) { return b1->getNumberOfStates() < b2->getNumberOfStates(); } );
-                
+
                 // Get and prepare the next splitter.
                 Block<BlockDataType>* splitter = splitterQueue.front();
                 splitterQueue.pop_front();

src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp

@@ -159,6 +159,8 @@ namespace storm {
         template<typename ModelType>
         void DeterministicModelBisimulationDecomposition<ModelType>::refinePredecessorBlocksOfSplitterStrong(std::list<Block<BlockDataType>*> const& predecessorBlocks, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
             for (auto block : predecessorBlocks) {
+                STORM_LOG_TRACE("Refining predecessor " << block->getId() << " of splitter");
+
                 // Depending on the actions we need to take, the block to refine changes, so we need to keep track of it.
                 Block<BlockDataType>* blockToRefineProbabilistically = block;
                 
@@ -428,6 +430,8 @@ namespace storm {
         
         template<typename ModelType>
         void DeterministicModelBisimulationDecomposition<ModelType>::refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue) {
+            STORM_LOG_TRACE("Refining partition based on splitter " << splitter.getId());
+            
             // The outline of the refinement is as follows.
             //
             // We iterate over all states of the splitter and determine for each predecessor the state the probability

src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp

@@ -74,9 +74,8 @@ namespace storm {
         
         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");
+            // TODO
         }
         
         template<typename ModelType>
@@ -129,17 +128,105 @@ namespace storm {
         
         template<typename ModelType>
         bool NondeterministicModelBisimulationDecomposition<ModelType>::checkQuotientDistributions() const {
-            
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (auto state = 0; state < this->model.getNumberOfStates(); ++state) {
+                for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                    storm::storage::Distribution<ValueType> distribution;
+                    for (auto const& element : this->model.getTransitionMatrix().getRow(choice)) {
+                        distribution.addProbability(this->partition.getBlock(element.getColumn()).getId(), element.getValue());
+                    }
+
+                    if (!distribution.equals(quotientDistributions[choice])) {
+                        std::cout << "the distributions for choice " << choice << " of state " << state << " do not match." << std::endl;
+                        std::cout << "is: " << quotientDistributions[choice] << " but should be " << distribution << std::endl;
+                        exit(-1);
+                    }
+
+                    bool less1 = distribution.less(quotientDistributions[choice], this->comparator);
+                    bool less2 = quotientDistributions[choice].less(distribution, this->comparator);
+                    
+                    if (distribution.equals(quotientDistributions[choice]) && (less1 || less2)) {
+                        std::cout << "mismatch of equality and less for " << std::endl;
+                        std::cout << quotientDistributions[choice] << " vs " << distribution << std::endl;
+                        exit(-1);
+                    }
+                }
+                
+                for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1] - 1; ++choice) {
+                    if (orderedQuotientDistributions[choice + 1]->less(*orderedQuotientDistributions[choice], this->comparator)) {
+                        std::cout << "choice " << (choice+1) << " is less than predecessor" << std::endl;
+                        std::cout << *orderedQuotientDistributions[choice] << " should be less than " << *orderedQuotientDistributions[choice + 1] << std::endl;
+                        exit(-1);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::printDistributions(uint_fast64_t state) const {
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (auto choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                std::cout << quotientDistributions[choice] << std::endl;
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkBlockStable(bisimulation::Block<BlockDataType> const& newBlock) const {
+            std::cout << "checking stability of new block " << newBlock.getId() << " of size " << newBlock.getNumberOfStates() << std::endl;
+            for (auto stateIt1 = this->partition.begin(newBlock), stateIte1 = this->partition.end(newBlock); stateIt1 != stateIte1; ++stateIt1) {
+                for (auto stateIt2 = this->partition.begin(newBlock), stateIte2 = this->partition.end(newBlock); stateIt2 != stateIte2; ++stateIt2) {
+                    bool less1 = quotientDistributionsLess(*stateIt1, *stateIt2);
+                    bool less2 = quotientDistributionsLess(*stateIt2, *stateIt1);
+                    if (less1 || less2) {
+                        std::cout << "the partition is not stable for the states " << *stateIt1 << " and " << *stateIt2 << std::endl;
+                        std::cout << "less1 " << less1 << " and less2 " << less2 << std::endl;
+                        
+                        std::cout << "distributions of state " << *stateIt1 << std::endl;
+                        this->printDistributions(*stateIt1);
+                        std::cout << "distributions of state " << *stateIt2 << std::endl;
+                        this->printDistributions(*stateIt2);
+                        exit(-1);
+                    }
+                }
+            }
+            return true;
+        }
+        
+        template<typename ModelType>
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::checkDistributionsDifferent(bisimulation::Block<BlockDataType> const& block, storm::storage::sparse::state_type end) const {
+            for (auto stateIt1 = this->partition.begin(block), stateIte1 = this->partition.end(block); stateIt1 != stateIte1; ++stateIt1) {
+                for (auto stateIt2 = this->partition.begin() + block.getEndIndex(), stateIte2 = this->partition.begin() + end; stateIt2 != stateIte2; ++stateIt2) {
+                    if (!quotientDistributionsLess(*stateIt1, *stateIt2)) {
+                        std::cout << "distributions are not less, even though they should be!" << std::endl;
+                        exit(-3);
+                    } else {
+                        std::cout << "less:" << std::endl;
+                        this->printDistributions(*stateIt1);
+                        std::cout << "and" << std::endl;
+                        this->printDistributions(*stateIt2);
+                    }
+                }
+            }
+            return true;
         }
         
         template<typename ModelType>
         bool NondeterministicModelBisimulationDecomposition<ModelType>::splitBlockAccordingToCurrentQuotientDistributions(Block<BlockDataType>& block, std::deque<Block<BlockDataType>*>& splitterQueue) {
+            std::list<Block<BlockDataType>*> newBlocks;
             bool split = this->partition.splitBlock(block,
                                                     [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
-                                                        return quotientDistributionsLess(state1, state2);
+                                                        bool result = quotientDistributionsLess(state1, state2);
+//                                                        std::cout << state1 << " is " << (!result ? "not" : "") << " smaller than " << state2 << std::endl;
+                                                        return result;
                                                     },
-                                                    [this, &block, &splitterQueue] (Block<BlockDataType>& newBlock) {
-                                                        updateQuotientDistributionsOfPredecessors(newBlock, block, splitterQueue);
+                                                    [this, &block, &splitterQueue, &newBlocks] (Block<BlockDataType>& newBlock) {
+                                                        newBlocks.push_back(&newBlock);
+                                                        
+//                                                        this->checkBlockStable(newBlock);
+//                                                        this->checkDistributionsDifferent(block, block.getEndIndex());
+//                                                        this->checkQuotientDistributions();
                                                     });
             
             // The quotient distributions of the predecessors of block do not need to be updated, since the probability
@@ -148,25 +235,26 @@ namespace storm {
 //            std::cout << "partition after split: " << std::endl;
 //            this->partition.print();
             
-            this->checkQuotientDistributions();
+            // defer updating the quotient distributions until *after* all splits, because
+            // it otherwise influences the subsequent splits!
+            for (auto el : newBlocks) {
+                this->updateQuotientDistributionsOfPredecessors(*el, block, splitterQueue);
+            }
+            
+//            this->checkQuotientDistributions();
             
             return split;
         }
         
         template<typename ModelType>
-        bool NondeterministicModelBisimulationDecomposition<ModelType>::quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) {
+        bool NondeterministicModelBisimulationDecomposition<ModelType>::quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) const {
             STORM_LOG_TRACE("Comparing the quotient distributions of state " << state1 << " and " << state2 << ".");
             std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
-            
-//            std::cout << "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;
-//            }
+
+//            std::cout << "distributions of state " << state1 << std::endl;
+//            this->printDistributions(state1);
+//            std::cout << "distributions of state " << state2 << std::endl;
+//            this->printDistributions(state2);
             
             auto firstIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1];
             auto firstIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1 + 1];

src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h

@@ -55,7 +55,7 @@ namespace storm {
             bool splitBlockAccordingToCurrentQuotientDistributions(bisimulation::Block<BlockDataType>& block, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
             
             // Retrieves whether the quotient distributions of state 1 are considered to be less than the ones of state 2.
-            bool quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2);
+            bool quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) const;
             
             // Updates the ordered list of quotient distribution for the given state.
             void updateOrderedQuotientDistributions(storm::storage::sparse::state_type state);
@@ -65,6 +65,9 @@ namespace storm {
             void updateQuotientDistributionsOfPredecessors(bisimulation::Block<BlockDataType> const& newBlock, bisimulation::Block<BlockDataType> const& oldBlock, std::deque<bisimulation::Block<BlockDataType>*>& splitterQueue);
             
             bool checkQuotientDistributions() const;
+            bool checkBlockStable(bisimulation::Block<BlockDataType> const& newBlock) const;
+            bool printDistributions(storm::storage::sparse::state_type state) const;
+            bool checkDistributionsDifferent(bisimulation::Block<BlockDataType> const& block, storm::storage::sparse::state_type end) const;
             
             // A mapping from choice indices to the state state that has this choice.
             std::vector<storm::storage::sparse::state_type> choiceToStateMapping;

src/storage/bisimulation/Partition.cpp

@@ -180,6 +180,7 @@ namespace storm {
             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);
                 }
@@ -210,6 +211,7 @@ namespace storm {
             template<typename DataType>
             std::pair<typename std::vector<std::unique_ptr<Block<DataType>>>::iterator, bool> Partition<DataType>::splitBlock(Block<DataType>& block, storm::storage::sparse::state_type position) {
                 STORM_LOG_THROW(position >= block.getBeginIndex() && position <= block.getEndIndex(), storm::exceptions::InvalidArgumentException, "Cannot split block at illegal position.");
+                STORM_LOG_TRACE("Splitting " << block.getId() << " at position " << position << " (begin was " << block.getBeginIndex() << ".");
 
                 // In case one of the resulting blocks would be empty, we simply return the current block and do not create
                 // a new one.
@@ -247,6 +249,7 @@ namespace storm {
                 std::vector<storm::storage::sparse::state_type>::const_iterator upperBound;
                 do {
                     upperBound = std::upper_bound(it, ite, *it, less);
+                    
                     if (upperBound != ite) {
                         wasSplit = true;
                         auto result = this->splitBlock(block, std::distance(this->states.cbegin(), upperBound));

src/utility/ConstantsComparator.cpp

@@ -101,7 +101,7 @@ namespace storm {
         }
         
         bool ConstantsComparator<double>::isLess(double const& value1, double const& value2) const {
-            return std::abs(value1 - value2) < precision;
+            return value1 < value2 - precision;
         }
         
         // Explicit instantiations.