extracting the bisimulation quotient for MDPs; tests for MDP bisimulation

Former-commit-id: 5613c653ba

src/builder/ExplicitPrismModelBuilder.h

@@ -72,8 +72,6 @@ namespace storm {
                 std::string stateInfo(uint_fast64_t state) const override {
                     return valuations[state].toString();
                 }
-
-
             };
             
             // A structure storing information about the used variables of the program.

src/storage/Distribution.cpp

@@ -32,11 +32,9 @@ 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;
                 }
             }

src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp

@@ -49,7 +49,20 @@ namespace storm {
         template<typename ModelType>
         void NondeterministicModelBisimulationDecomposition<ModelType>::initializeQuotientDistributions() {
             std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
-            for (auto choice = 0; choice < nondeterministicChoiceIndices.back(); ++choice) {
+            
+            for (auto const& block : this->partition.getBlocks()) {
+                if (block->data().absorbing()) {
+                    // If the block is marked as absorbing, we need to create the corresponding distributions.
+                    for (auto stateIt = this->partition.begin(*block), stateIte = this->partition.end(*block); stateIt != stateIte; ++stateIt) {
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[*stateIt]; choice < nondeterministicChoiceIndices[*stateIt + 1]; ++choice) {
+                            this->quotientDistributions[choice].addProbability(block->getId(), storm::utility::one<ValueType>());
+                            orderedQuotientDistributions[choice] = &this->quotientDistributions[choice];
+                        }
+                    }
+                } else {
+                    // Otherwise, we compute the probabilities from the transition matrix.
+                    for (auto stateIt = this->partition.begin(*block), stateIte = this->partition.end(*block); stateIt != stateIte; ++stateIt) {
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[*stateIt]; choice < nondeterministicChoiceIndices[*stateIt + 1]; ++choice) {
                             for (auto entry : this->model.getTransitionMatrix().getRow(choice)) {
                                 if (!this->comparator.isZero(entry.getValue())) {
                                     this->quotientDistributions[choice].addProbability(this->partition.getBlock(entry.getColumn()).getId(), entry.getValue());
@@ -57,6 +70,9 @@ namespace storm {
                             }
                             orderedQuotientDistributions[choice] = &this->quotientDistributions[choice];
                         }
+                    }
+                }
+            }
             
             for (auto state = 0; state < this->model.getNumberOfStates(); ++state) {
                 updateOrderedQuotientDistributions(state);
@@ -74,8 +90,106 @@ namespace storm {
         
         template<typename ModelType>
         void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
-            STORM_LOG_ASSERT(false, "Not yet implemented");
-            // TODO
+            // In order to create the quotient model, we need to construct
+            // (a) the new transition matrix,
+            // (b) the new labeling,
+            // (c) the new reward structures.
+            
+            // Prepare a matrix builder for (a).
+            storm::storage::SparseMatrixBuilder<ValueType> builder(0, this->size(), 0, false, true, this->size());
+            
+            // Prepare the new state labeling for (b).
+            storm::models::sparse::StateLabeling newLabeling(this->size());
+            std::set<std::string> atomicPropositionsSet = this->options.respectedAtomicPropositions.get();
+            atomicPropositionsSet.insert("init");
+            std::vector<std::string> atomicPropositions = std::vector<std::string>(atomicPropositionsSet.begin(), atomicPropositionsSet.end());
+            for (auto const& ap : atomicPropositions) {
+                newLabeling.addLabel(ap);
+            }
+            
+            // If the model had state rewards, we need to build the state rewards for the quotient as well.
+            boost::optional<std::vector<ValueType>> stateRewards;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                stateRewards = std::vector<ValueType>(this->blocks.size());
+            }
+            
+            // Now build (a) and (b) by traversing all blocks.
+            uint_fast64_t currentRow = 0;
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices();
+            for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) {
+                auto const& block = this->blocks[blockIndex];
+                
+                // Open new row group for the new meta state.
+                builder.newRowGroup(currentRow);
+                
+                // Pick one representative state. For strong bisimulation it doesn't matter which state it is, because
+                // they all behave equally.
+                storm::storage::sparse::state_type representativeState = *block.begin();
+                Block<BlockDataType> const& oldBlock = this->partition.getBlock(representativeState);
+                
+                // If the block is absorbing, we simply add a self-loop.
+                if (oldBlock.data().absorbing()) {
+                    builder.addNextValue(currentRow, blockIndex, storm::utility::one<ValueType>());
+                    ++currentRow;
+                    
+                    // If the block has a special representative state, we retrieve it now.
+                    if (oldBlock.data().hasRepresentativeState()) {
+                        representativeState = oldBlock.data().representativeState();
+                    }
+                    
+                    // Add all of the selected atomic propositions that hold in the representative state to the state
+                    // representing the block.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                } else {
+                    // Add the outgoing choices of the block.
+                    for (uint_fast64_t choice = nondeterministicChoiceIndices[representativeState]; choice < nondeterministicChoiceIndices[representativeState + 1]; ++choice) {
+                        // If the choice is the same as the last one, we do not need to add it.
+                        if (choice > nondeterministicChoiceIndices[representativeState] && quotientDistributions[choice - 1].equals(quotientDistributions[choice], this->comparator)) {
+                            continue;
+                        }
+                        
+                        for (auto entry : quotientDistributions[choice]) {
+                            builder.addNextValue(currentRow, entry.first, entry.second);
+                        }
+                        ++currentRow;
+                    }
+                    
+                    // Otherwise add all atomic propositions to the equivalence class that the representative state
+                    // satisfies.
+                    for (auto const& ap : atomicPropositions) {
+                        if (this->model.getStateLabeling().getStateHasLabel(ap, representativeState)) {
+                            newLabeling.addLabelToState(ap, blockIndex);
+                        }
+                    }
+                }
+                
+                // If the model has state rewards, we simply copy the state reward of the representative state, because
+                // all states in a block are guaranteed to have the same state reward.
+                if (this->options.keepRewards && this->model.hasRewardModel()) {
+                    typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                    stateRewards.get()[blockIndex] = nameRewardModelPair->second.getStateRewardVector()[representativeState];
+                }
+            }
+            
+            // Now check which of the blocks of the partition contain at least one initial state.
+            for (auto initialState : this->model.getInitialStates()) {
+                Block<BlockDataType> const& initialBlock = this->partition.getBlock(initialState);
+                newLabeling.addLabelToState("init", initialBlock.getId());
+            }
+            
+            // Construct the reward model mapping.
+            std::unordered_map<std::string, typename ModelType::RewardModelType> rewardModels;
+            if (this->options.keepRewards && this->model.hasRewardModel()) {
+                typename std::unordered_map<std::string, typename ModelType::RewardModelType>::const_iterator nameRewardModelPair = this->model.getUniqueRewardModel();
+                rewardModels.insert(std::make_pair(nameRewardModelPair->first, typename ModelType::RewardModelType(stateRewards)));
+            }
+            
+            // Finally construct the quotient model.
+            this->quotient = std::shared_ptr<ModelType>(new ModelType(builder.build(), std::move(newLabeling), std::move(rewardModels)));
         }
         
         template<typename ModelType>
@@ -98,6 +212,11 @@ namespace storm {
                     storm::storage::sparse::state_type predecessorState = choiceToStateMapping[predecessorChoice];
                     Block<BlockDataType>& predecessorBlock = this->partition.getBlock(predecessorState);
                     
+                    // If the predecessor block is marked as absorbing, we do not need to update anything.
+                    if (predecessorBlock.data().absorbing()) {
+                        continue;
+                    }
+                    
                     // If the predecessor block is not marked as to-refined, we do so now.
                     if (!predecessorBlock.data().splitter()) {
                         predecessorBlock.data().setSplitter();

test/functional/storage/NondeterministicModelBisimulationDecompositionTest.cpp

@@ -0,0 +1,58 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/parser/PrismParser.h"
+#include "src/parser/FormulaParser.h"
+
+#include "src/builder/ExplicitPrismModelBuilder.h"
+
+#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
+#include "src/models/sparse/Mdp.h"
+#include "src/models/sparse/StandardRewardModel.h"
+
+TEST(NondeterministicModelBisimulationDecomposition, TwoDice) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/two_dice.nm");
+
+    // Build the die model without its reward model.
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    std::shared_ptr<storm::models::sparse::Mdp<double>> mdp = model->as<storm::models::sparse::Mdp<double>>();
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim(*mdp);
+    ASSERT_NO_THROW(bisim.computeBisimulationDecomposition());
+    std::shared_ptr<storm::models::sparse::Model<double>> result;
+    ASSERT_NO_THROW(result = bisim.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(77ul, result->getNumberOfStates());
+    EXPECT_EQ(183ul, result->getNumberOfTransitions());
+    EXPECT_EQ(97ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options;
+    options.respectedAtomicPropositions = std::set<std::string>({"two"});
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim2(*mdp, options);
+    ASSERT_NO_THROW(bisim2.computeBisimulationDecomposition());
+    ASSERT_NO_THROW(result = bisim2.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(11ul, result->getNumberOfStates());
+    EXPECT_EQ(26ul, result->getNumberOfTransitions());
+    EXPECT_EQ(14ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    std::shared_ptr<storm::logic::Formula> formula = formulaParser.parseSingleFormulaFromString("Pmin=? [F \"two\"]");
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options2(*mdp, *formula);
+    
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim3(*mdp, options2);
+    ASSERT_NO_THROW(bisim3.computeBisimulationDecomposition());
+    ASSERT_NO_THROW(result = bisim3.getQuotient());
+    
+    EXPECT_EQ(storm::models::ModelType::Mdp, result->getType());
+    EXPECT_EQ(11ul, result->getNumberOfStates());
+    EXPECT_EQ(26ul, result->getNumberOfTransitions());
+    EXPECT_EQ(14ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
+}