more work on symbolic bisimulation

resources/3rdparty/CMakeLists.txt

@@ -383,7 +383,7 @@ ExternalProject_Add(
         LOG_CONFIGURE ON
         LOG_BUILD ON
         BUILD_BYPRODUCTS ${STORM_3RDPARTY_BINARY_DIR}/sylvan/src/libsylvan${STATIC_EXT}
-        BUILD_ALWAYS
+        BUILD_ALWAYS 1
 )
 ExternalProject_Get_Property(sylvan source_dir)
 ExternalProject_Get_Property(sylvan binary_dir)

src/storm/storage/dd/BisimulationDecomposition.cpp

@@ -0,0 +1,84 @@
+#include "storm/storage/dd/BisimulationDecomposition.h"
+
+#include "storm/storage/dd/bisimulation/SignatureRefiner.h"
+
+#include "storm/utility/macros.h"
+#include "storm/exceptions/InvalidOperationException.h"
+
+#include <sylvan_table.h>
+
+extern llmsset_t nodes;
+
+namespace storm {
+    namespace dd {
+        
+        using namespace bisimulation;
+        
+        template <storm::dd::DdType DdType, typename ValueType>
+        BisimulationDecomposition<DdType, ValueType>::BisimulationDecomposition(storm::models::symbolic::Model<DdType, ValueType> const& model) : status(Status::Initialized), model(model), currentPartition(bisimulation::Partition<DdType, ValueType>::create(model)) {
+            // Intentionally left empty.
+        }
+        
+        template <storm::dd::DdType DdType, typename ValueType>
+        BisimulationDecomposition<DdType, ValueType>::BisimulationDecomposition(storm::models::symbolic::Model<DdType, ValueType> const& model, Partition<DdType, ValueType> const& initialPartition) : model(model), currentPartition(initialPartition) {
+            // Intentionally left empty.
+        }
+        
+        template <storm::dd::DdType DdType, typename ValueType>
+        void BisimulationDecomposition<DdType, ValueType>::compute() {
+//            LACE_ME;
+            
+            auto partitionRefinementStart = std::chrono::high_resolution_clock::now();
+            this->status = Status::InComputation;
+            
+            STORM_LOG_TRACE("Initial partition has " << currentPartition.getNumberOfBlocks() << " blocks.");
+#ifndef NDEBUG
+            STORM_LOG_TRACE("Initial partition ADD has " << currentPartition.getPartitionAdd().getNodeCount() << " nodes.");
+#endif
+            
+            SignatureRefiner<DdType, ValueType> refiner(model.getManager(), currentPartition.getBlockVariable(), model.getRowVariables());
+            bool done = false;
+            uint64_t iterations = 0;
+            while (!done) {
+//                currentPartition.getPartitionAdd().exportToDot("part" + std::to_string(iterations) + ".dot");
+                Signature<DdType, ValueType> signature(model.getTransitionMatrix().multiplyMatrix(currentPartition.getPartitionAdd(), model.getColumnVariables()));
+//                signature.getSignatureAdd().exportToDot("sig" + std::to_string(iterations) + ".dot");
+#ifndef NDEBUG
+                STORM_LOG_TRACE("Computed signature ADD with " << signature.getSignatureAdd().getNodeCount() << " nodes.");
+#endif
+                
+                Partition<DdType, ValueType> newPartition = refiner.refine(currentPartition, signature);
+
+                STORM_LOG_TRACE("New partition has " << newPartition.getNumberOfBlocks() << " blocks.");
+#ifndef NDEBUG
+                STORM_LOG_TRACE("Computed new partition ADD with " << newPartition.getPartitionAdd().getNodeCount() << " nodes.");
+#endif
+//                STORM_LOG_TRACE("Current #nodes in table " << llmsset_count_marked(nodes) << " of " << llmsset_get_size(nodes) << " BDD nodes.");
+                
+                if (currentPartition == newPartition) {
+                    done = true;
+                } else {
+                    currentPartition = newPartition;
+                }
+                ++iterations;
+            }
+
+            this->status = Status::FixedPoint;
+            auto partitionRefinementEnd = std::chrono::high_resolution_clock::now();
+            STORM_LOG_DEBUG("Partition refinement completed in " << std::chrono::duration_cast<std::chrono::milliseconds>(partitionRefinementEnd - partitionRefinementStart).count() << "ms (" << iterations << " iterations).");
+        }
+        
+        template <storm::dd::DdType DdType, typename ValueType>
+        std::shared_ptr<storm::models::symbolic::Model<DdType, ValueType>> BisimulationDecomposition<DdType, ValueType>::getQuotient() const {
+            STORM_LOG_THROW(this->status == Status::FixedPoint, storm::exceptions::InvalidOperationException, "Cannot extract quotient, because bisimulation decomposition was not completed.");
+            return nullptr;
+        }
+        
+        template class BisimulationDecomposition<storm::dd::DdType::CUDD, double>;
+
+        template class BisimulationDecomposition<storm::dd::DdType::Sylvan, double>;
+        template class BisimulationDecomposition<storm::dd::DdType::Sylvan, storm::RationalNumber>;
+        template class BisimulationDecomposition<storm::dd::DdType::Sylvan, storm::RationalFunction>;
+        
+    }
+}

src/storm/storage/dd/BisimulationDecomposition.h

@@ -0,0 +1,45 @@
+#pragma once
+
+#include "storm/storage/dd/DdType.h"
+
+#include "storm/models/symbolic/Model.h"
+
+#include "storm/storage/dd/bisimulation/Partition.h"
+
+namespace storm {
+    namespace dd {
+        
+        template <storm::dd::DdType DdType, typename ValueType>
+        class BisimulationDecomposition {
+        public:
+            enum class Status {
+                Initialized, InComputation, FixedPoint
+            };
+            
+            BisimulationDecomposition(storm::models::symbolic::Model<DdType, ValueType> const& model);
+
+            BisimulationDecomposition(storm::models::symbolic::Model<DdType, ValueType> const& model, bisimulation::Partition<DdType, ValueType> const& initialPartition);
+            
+            /*!
+             * Computes the decomposition.
+             */
+            void compute();
+            
+            /*!
+             * Retrieves the quotient model after the bisimulation decomposition was computed.
+             */
+            std::shared_ptr<storm::models::symbolic::Model<DdType, ValueType>> getQuotient() const;
+            
+        private:
+            // The status of the computation.
+            Status status;
+            
+            // The model for which to compute the bisimulation decomposition.
+            storm::models::symbolic::Model<DdType, ValueType> const& model;
+            
+            // The current partition in the partition refinement process. Initially set to the initial partition.
+            bisimulation::Partition<DdType, ValueType> currentPartition;
+        };
+        
+    }
+}

src/storm/storage/dd/DdManager.cpp

@@ -56,7 +56,7 @@ namespace storm {
         }
         
         template<DdType LibraryType>
-        Bdd<LibraryType> DdManager<LibraryType>::getEncoding(storm::expressions::Variable const& variable, int_fast64_t value) const {
+        Bdd<LibraryType> DdManager<LibraryType>::getEncoding(storm::expressions::Variable const& variable, int_fast64_t value, bool mostSignificantBitAtTop) const {
             DdMetaVariable<LibraryType> const& metaVariable = this->getMetaVariable(variable);
             
             STORM_LOG_THROW(value >= metaVariable.getLow() && value <= metaVariable.getHigh(), storm::exceptions::InvalidArgumentException, "Illegal value " << value << " for meta variable '" << variable.getName() << "'.");
@@ -67,17 +67,35 @@ namespace storm {
             std::vector<Bdd<LibraryType>> const& ddVariables = metaVariable.getDdVariables();
             
             Bdd<LibraryType> result;
-            if (value & (1ull << (ddVariables.size() - 1))) {
-                result = ddVariables[0];
+            if (mostSignificantBitAtTop) {
+                if (value & (1ull << (ddVariables.size() - 1))) {
+                    result = ddVariables[0];
+                } else {
+                    result = !ddVariables[0];
+                }
+                
+                for (std::size_t i = 1; i < ddVariables.size(); ++i) {
+                    if (value & (1ull << (ddVariables.size() - i - 1))) {
+                        result &= ddVariables[i];
+                    } else {
+                        result &= !ddVariables[i];
+                    }
+                }
             } else {
-                result = !ddVariables[0];
-            }
-            
-            for (std::size_t i = 1; i < ddVariables.size(); ++i) {
-                if (value & (1ull << (ddVariables.size() - i - 1))) {
-                    result &= ddVariables[i];
+                if (value & 1ull) {
+                    result = ddVariables[0];
                 } else {
-                    result &= !ddVariables[i];
+                    result = !ddVariables[0];
+                }
+                value >>= 1;
+                
+                for (std::size_t i = 1; i < ddVariables.size(); ++i) {
+                    if (value & 1ull) {
+                        result &= ddVariables[i];
+                    } else {
+                        result &= !ddVariables[i];
+                    }
+                    value >>= 1;
                 }
             }
             
@@ -164,6 +182,8 @@ namespace storm {
                 ++numberOfBits;
             }
             
+            STORM_LOG_TRACE("Creating meta variable with " << numberOfBits << " bit(s) and " << numberOfLayers << " layer(s).");
+            
             std::stringstream tmp1;
             std::vector<storm::expressions::Variable> result;
             for (uint64 layer = 0; layer < numberOfLayers; ++layer) {

src/storm/storage/dd/DdManager.h

@@ -94,10 +94,12 @@ namespace storm {
              *
              * @param variable The expression variable associated with the meta variable.
              * @param value The value the meta variable is supposed to have.
+             * @param mostSignificantBitAtTop A flag indicating whether the most significant bit of the value is to be
+             * encoded with the topmost variable or the bottommost.
              * @return The DD representing the function that maps all inputs which have the given meta variable equal
              * to the given value one.
              */
-            Bdd<LibraryType> getEncoding(storm::expressions::Variable const& variable, int_fast64_t value) const;
+            Bdd<LibraryType> getEncoding(storm::expressions::Variable const& variable, int_fast64_t value, bool mostSignificantBitAtTop = true) const;
             
             /*!
              * Retrieves the BDD representing the range of the meta variable, i.e., a function that maps all legal values

src/storm/storage/dd/bisimulation/Partition.cpp

@@ -90,10 +90,13 @@ namespace storm {
                 // Enumerate all realizable blocks.
                 enumerateBlocksRec<DdType>(stateSets, model.getReachableStates(), 0, blockVariable, [&manager, &partitionAdd, &blockVariable, &blockCount](storm::dd::Bdd<DdType> const& stateSet) {
                     stateSet.template toAdd<ValueType>().exportToDot("states_" + std::to_string(blockCount) + ".dot");
-                    partitionAdd += (stateSet && manager.getEncoding(blockVariable, blockCount)).template toAdd<ValueType>();
+                    partitionAdd += (stateSet && manager.getEncoding(blockVariable, blockCount, false)).template toAdd<ValueType>();
                     blockCount++;
                 } );
                 
+                // Move the partition over to the primed variables.
+                partitionAdd = partitionAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                
                 return std::make_pair(partitionAdd, blockCount);
             }
             

src/storm/storage/dd/bisimulation/SignatureRefiner.cpp

@@ -75,7 +75,7 @@ namespace storm {
                     }
                     
                     // Determine the levels in the DDs.
-                    uint64_t partitionVariable = Cudd_NodeReadIndex(partitionNode);
+                    uint64_t partitionVariable = Cudd_NodeReadIndex(partitionNode) - 1;
                     uint64_t signatureVariable = Cudd_NodeReadIndex(signatureNode);
                     uint64_t partitionLevel = Cudd_ReadPerm(ddman, partitionVariable);
                     uint64_t signatureLevel = Cudd_ReadPerm(ddman, signatureVariable);
@@ -111,7 +111,7 @@ namespace storm {
                         }
                         
                         // Get the node to connect the subresults.
-                        DdNode* var = Cudd_addIthVar(ddman, topVariable);
+                        DdNode* var = Cudd_addIthVar(ddman, topVariable + 1);
                         Cudd_Ref(var);
                         DdNode* result = Cudd_addIte(ddman, var, thenResult, elseResult);
                         Cudd_Ref(result);
@@ -144,7 +144,7 @@ namespace storm {
                         } else {
                             DdNode* result;
                             {
-                                storm::dd::Add<storm::dd::DdType::CUDD, ValueType> blockEncoding = manager.getEncoding(blockVariable, nextFreeBlockIndex).template toAdd<ValueType>();
+                                storm::dd::Add<storm::dd::DdType::CUDD, ValueType> blockEncoding = manager.getEncoding(blockVariable, nextFreeBlockIndex, false).template toAdd<ValueType>();
                                 ++nextFreeBlockIndex;
                                 result = blockEncoding.getInternalAdd().getCuddDdNode();
                                 Cudd_Ref(result);
@@ -175,8 +175,7 @@ namespace storm {
                 // The cache used to identify which old block numbers have already been reused.
                 spp::sparse_hash_map<DdNode const*, bool> reuseBlocksCache;
             };
-            
-            
+                        
             template<typename ValueType>
             class InternalSignatureRefiner<storm::dd::DdType::Sylvan, ValueType> {
             public:
@@ -219,7 +218,7 @@ namespace storm {
 
                     // Determine levels in the DDs.
                     BDDVAR signatureVariable = mtbdd_isleaf(signatureNode) ? 0xffffffff : sylvan_var(signatureNode);
-                    BDDVAR partitionVariable = sylvan_var(partitionNode);
+                    BDDVAR partitionVariable = mtbdd_isleaf(signatureNode) ? 0xffffffff : sylvan_var(partitionNode) - 1;
                     BDDVAR topVariable = std::min(signatureVariable, partitionVariable);
 
                     // Check whether the top variable is still within the state encoding.
@@ -229,25 +228,37 @@ namespace storm {
                         MTBDD elseResult;
                         if (partitionVariable < signatureVariable) {
                             thenResult = refine(sylvan_high(partitionNode), signatureNode);
+                            sylvan_protect(&thenResult);
                             elseResult = refine(sylvan_low(partitionNode), signatureNode);
+                            sylvan_protect(&elseResult);
                         } else if (partitionVariable > signatureVariable) {
                             thenResult = refine(partitionNode, sylvan_high(signatureNode));
+                            sylvan_protect(&thenResult);
                             elseResult = refine(partitionNode, sylvan_low(signatureNode));
+                            sylvan_protect(&elseResult);
                         } else {
                             thenResult = refine(sylvan_high(partitionNode), sylvan_high(signatureNode));
+                            sylvan_protect(&thenResult);
                             elseResult = refine(sylvan_low(partitionNode), sylvan_low(signatureNode));
+                            sylvan_protect(&elseResult);
                         }
                         
                         if (thenResult == elseResult) {
+                            sylvan_unprotect(&thenResult);
+                            sylvan_unprotect(&elseResult);
                             return thenResult;
                         }
                         
                         // Get the node to connect the subresults.
-                        MTBDD result = mtbdd_ite(sylvan_ithvar(topVariable), thenResult, elseResult);
+                        MTBDD result = sylvan_makenode(topVariable + 1, elseResult, thenResult);// mtbdd_ite(sylvan_ithvar(topVariable), thenResult, elseResult);
+                        sylvan_protect(&result);
+                        sylvan_unprotect(&thenResult);
+                        sylvan_unprotect(&elseResult);
                         
                         // Store the result in the cache.
                         signatureCache[std::make_pair(signatureNode, partitionNode)] = result;
-                        
+
+                        sylvan_unprotect(&result);
                         return result;
                     } else {
                         
@@ -269,11 +280,13 @@ namespace storm {
                         } else {
                             MTBDD result;
                             {
-                                storm::dd::Add<storm::dd::DdType::Sylvan, ValueType> blockEncoding = manager.getEncoding(blockVariable, nextFreeBlockIndex).template toAdd<ValueType>();
+                                storm::dd::Add<storm::dd::DdType::Sylvan, ValueType> blockEncoding = manager.getEncoding(blockVariable, nextFreeBlockIndex, false).template toAdd<ValueType>();
                                 ++nextFreeBlockIndex;
                                 result = blockEncoding.getInternalAdd().getSylvanMtbdd().GetMTBDD();
+                                sylvan_protect(&result);
                             }
                             signatureCache[std::make_pair(signatureNode, partitionNode)] = result;
+                            sylvan_unprotect(&result);
                             return result;
                         }
                     }

src/storm/storage/dd/sylvan/InternalSylvanDdManager.cpp

@@ -36,9 +36,9 @@ namespace storm {
             if (numberOfInstances == 0) {
                 storm::settings::modules::SylvanSettings const& settings = storm::settings::getModule<storm::settings::modules::SylvanSettings>();
                 if (settings.isNumberOfThreadsSet()) {
-                    lace_init(settings.getNumberOfThreads(), 1000000);
+                    lace_init(settings.getNumberOfThreads(), 1024*1024*16);
                 } else {
-                    lace_init(0, 1000000);
+                    lace_init(0, 1024*1024*16);
                 }
                 lace_startup(0, 0, 0);
                 
@@ -50,8 +50,11 @@ namespace storm {
                 
                 STORM_LOG_THROW(powerOfTwo >= 16, storm::exceptions::InvalidSettingsException, "Too little memory assigned to sylvan.");
                 
-                STORM_LOG_TRACE("Assigning " << (1ull << (powerOfTwo - 1)) << " slots to both sylvan's unique table and its cache.");
-                sylvan::Sylvan::initPackage(1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1));
+                if ((((1ull << powerOfTwo) + (1ull << (powerOfTwo - 1)))) < totalNodesToStore) {
+                    sylvan::Sylvan::initPackage(1ull << powerOfTwo, 1ull << powerOfTwo, 1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1));
+                } else {
+                    sylvan::Sylvan::initPackage(1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1), 1ull << (powerOfTwo - 1));
+                }
                 sylvan::Sylvan::setGranularity(3);
                 sylvan::Sylvan::initBdd();
                 sylvan::Sylvan::initMtbdd();