diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 79b53b408..49d0cddad 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -1603,7 +1603,6 @@ namespace storm { for (index_type row = 0; row < this->rowCount; ++row) { auto rowSum = getRowSum(row); if (!comparator.isOne(rowSum)) { - std::cout << "row sum of row " << row << " is " << rowSum << std::endl; return false; } } diff --git a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp b/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp index 2a52fde56..fc6e289b1 100644 --- a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp +++ b/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp @@ -32,6 +32,182 @@ namespace storm { namespace dd { namespace bisimulation { + template + class InternalRepresentativeComputer; + + template + class InternalRepresentativeComputerBase { + public: + InternalRepresentativeComputerBase(Partition const& partition, std::set const& rowVariables, std::set const& columnVariables) : partition(partition), rowVariables(rowVariables), columnVariables(columnVariables) { + if (partition.storedAsAdd()) { + ddManager = &partition.asAdd().getDdManager(); + } else { + ddManager = &partition.asBdd().getDdManager(); + } + internalDdManager = &ddManager->getInternalDdManager(); + + // Create state variables cube. + this->columnVariablesCube = ddManager->getBddOne(); + for (auto const& var : columnVariables) { + auto const& metaVariable = ddManager->getMetaVariable(var); + this->columnVariablesCube &= metaVariable.getCube(); + } + } + + protected: + storm::dd::DdManager const* ddManager; + storm::dd::InternalDdManager const* internalDdManager; + + Partition const& partition; + std::set const& rowVariables; + std::set const& columnVariables; + storm::dd::Bdd columnVariablesCube; + }; + + template + class InternalRepresentativeComputer : public InternalRepresentativeComputerBase { + public: + InternalRepresentativeComputer(Partition const& partition, std::set const& rowVariables, std::set const& columnVariables) : InternalRepresentativeComputerBase(partition, rowVariables, columnVariables) { + this->ddman = this->internalDdManager->getCuddManager().getManager(); + } + + storm::dd::Bdd getRepresentatives() { + return storm::dd::Bdd(*this->ddManager, storm::dd::InternalBdd(this->internalDdManager, cudd::BDD(this->internalDdManager->getCuddManager(), this->getRepresentativesRec(this->partition.asAdd().getInternalAdd().getCuddDdNode(), this->columnVariablesCube.getInternalBdd().getCuddDdNode()))), this->rowVariables); + } + + private: + DdNodePtr getRepresentativesRec(DdNodePtr partitionNode, DdNodePtr stateVariablesCube) { + if (partitionNode == Cudd_ReadZero(ddman)) { + return Cudd_ReadLogicZero(ddman); + } + + // If we visited the node before, there is no block that we still need to cover. + if (visitedNodes.find(partitionNode) != visitedNodes.end()) { + return Cudd_ReadLogicZero(ddman); + } + + // If we hit a block variable and have not yet terminated the DFS earlier, it means we have a new representative. + if (Cudd_IsConstant(stateVariablesCube)) { + visitedNodes.emplace(partitionNode, true); + return Cudd_ReadOne(ddman); + } else { + bool skipped = false; + DdNodePtr elsePartitionNode; + DdNodePtr thenPartitionNode; + if (Cudd_NodeReadIndex(partitionNode) == Cudd_NodeReadIndex(stateVariablesCube)) { + elsePartitionNode = Cudd_E(partitionNode); + thenPartitionNode = Cudd_T(partitionNode); + } else { + elsePartitionNode = thenPartitionNode = partitionNode; + skipped = true; + } + + if (!skipped) { + visitedNodes.emplace(partitionNode, true); + } + + // Otherwise, recursively proceed with DFS. + DdNodePtr elseResult = getRepresentativesRec(elsePartitionNode, Cudd_T(stateVariablesCube)); + Cudd_Ref(elseResult); + + DdNodePtr thenResult = nullptr; + if (!skipped) { + thenResult = getRepresentativesRec(thenPartitionNode, Cudd_T(stateVariablesCube)); + Cudd_Ref(thenResult); + + if (thenResult == elseResult) { + Cudd_Deref(elseResult); + Cudd_Deref(thenResult); + return elseResult; + } else { + bool complement = Cudd_IsComplement(thenResult); + auto result = cuddUniqueInter(ddman, Cudd_NodeReadIndex(stateVariablesCube) - 1, Cudd_Regular(thenResult), complement ? Cudd_Not(elseResult) : elseResult); + Cudd_Deref(elseResult); + Cudd_Deref(thenResult); + return complement ? Cudd_Not(result) : result; + } + } else { + auto result = Cudd_Not(cuddUniqueInter(ddman, Cudd_NodeReadIndex(stateVariablesCube) - 1, Cudd_ReadOne(ddman), Cudd_Not(elseResult))); + Cudd_Deref(elseResult); + return result; + } + } + } + + ::DdManager* ddman; + spp::sparse_hash_map visitedNodes; + }; + + template + class InternalRepresentativeComputer : public InternalRepresentativeComputerBase { + public: + InternalRepresentativeComputer(Partition const& partition, std::set const& rowVariables, std::set const& columnVariables) : InternalRepresentativeComputerBase(partition, rowVariables, columnVariables) { + // Intentionally left empty. + } + + storm::dd::Bdd getRepresentatives() { + return storm::dd::Bdd(*this->ddManager, storm::dd::InternalBdd(this->internalDdManager, sylvan::Bdd(this->getRepresentativesRec(this->partition.asBdd().getInternalBdd().getSylvanBdd().GetBDD(), this->columnVariablesCube.getInternalBdd().getSylvanBdd().GetBDD()))), this->rowVariables); + } + + private: + BDD getRepresentativesRec(BDD partitionNode, BDD stateVariablesCube) { + if (partitionNode == sylvan_false) { + return sylvan_false; + } + + // If we visited the node before, there is no block that we still need to cover. + if (visitedNodes.find(partitionNode) != visitedNodes.end()) { + return sylvan_false; + } + + // If we hit a block variable and have not yet terminated the DFS earlier, it means we have a new representative. + if (sylvan_isconst(stateVariablesCube)) { + visitedNodes.emplace(partitionNode, true); + return sylvan_true; + } else { + bool skipped = false; + BDD elsePartitionNode; + BDD thenPartitionNode; + if (sylvan_var(partitionNode) == sylvan_var(stateVariablesCube)) { + elsePartitionNode = sylvan_low(partitionNode); + thenPartitionNode = sylvan_high(partitionNode); + } else { + elsePartitionNode = thenPartitionNode = partitionNode; + skipped = true; + } + + if (!skipped) { + visitedNodes.emplace(partitionNode, true); + } + + // Otherwise, recursively proceed with DFS. + BDD elseResult = getRepresentativesRec(elsePartitionNode, sylvan_high(stateVariablesCube)); + mtbdd_refs_push(elseResult); + + BDD thenResult; + if (!skipped) { + thenResult = getRepresentativesRec(thenPartitionNode, sylvan_high(stateVariablesCube)); + mtbdd_refs_push(thenResult); + + if (thenResult == elseResult) { + mtbdd_refs_pop(2); + return elseResult; + } else { + auto result = sylvan_makenode(sylvan_var(stateVariablesCube) - 1, elseResult, thenResult); + mtbdd_refs_pop(2); + return result; + } + } else { + auto result = sylvan_makenode(sylvan_var(stateVariablesCube) - 1, elseResult, sylvan_false); + mtbdd_refs_pop(1); + return result; + } + } + } + + spp::sparse_hash_map visitedNodes; + }; + template class InternalSparseQuotientExtractor; @@ -651,14 +827,23 @@ namespace storm { InternalSparseQuotientExtractor sparseExtractor(model.getManager(), model.getRowVariables(), model.getNondeterminismVariables()); auto states = partition.getStates().swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Bdd partitionAsBdd = partition.storedAsAdd() ? partition.asAdd().toBdd() : partition.asBdd(); + partitionAsBdd = partitionAsBdd.renameVariables(model.getColumnVariables(), model.getRowVariables()); + partitionAsBdd.template toAdd().exportToDot("partition.dot"); + auto start = std::chrono::high_resolution_clock::now(); - storm::storage::SparseMatrix quotientTransitionMatrix = sparseExtractor.extractTransitionMatrix(model.getTransitionMatrix(), partition); + // FIXME: Use partition as BDD in representative computation. + auto representatives = InternalRepresentativeComputer(partition, model.getRowVariables(), model.getColumnVariables()).getRepresentatives(); + representatives.template toAdd().exportToDot("repr.dot"); + (model.getTransitionMatrix() * representatives.template toAdd()).exportToDot("extract.dot"); + storm::storage::SparseMatrix quotientTransitionMatrix = sparseExtractor.extractTransitionMatrix(model.getTransitionMatrix() * representatives.template toAdd(), partition); auto end = std::chrono::high_resolution_clock::now(); STORM_LOG_TRACE("Quotient transition matrix extracted in " << std::chrono::duration_cast(end - start).count() << "ms."); start = std::chrono::high_resolution_clock::now(); + storm::dd::Odd odd = representatives.createOdd(); storm::models::sparse::StateLabeling quotientStateLabeling(partition.getNumberOfBlocks()); - quotientStateLabeling.addLabel("init", sparseExtractor.extractStates(model.getInitialStates(), partition)); + quotientStateLabeling.addLabel("init", ((model.getInitialStates() && partitionAsBdd).existsAbstract(model.getRowVariables()) && partitionAsBdd && representatives).existsAbstract({partition.getBlockVariable()}).toVector(odd)); quotientStateLabeling.addLabel("deadlock", sparseExtractor.extractStates(model.getDeadlockStates(), partition)); for (auto const& label : preservationInformation.getLabels()) { @@ -743,15 +928,15 @@ namespace storm { auto end = std::chrono::high_resolution_clock::now(); STORM_LOG_TRACE("Quotient labels extracted in " << std::chrono::duration_cast(end - start).count() << "ms."); - storm::dd::Add partitionAsAdd = partitionAsBdd.template toAdd(); start = std::chrono::high_resolution_clock::now(); - storm::dd::Add quotientTransitionMatrix = model.getTransitionMatrix().multiplyMatrix(partitionAsAdd.renameVariables(blockVariableSet, blockPrimeVariableSet), model.getColumnVariables()); + storm::dd::Add quotientTransitionMatrix = model.getTransitionMatrix().multiplyMatrix(partitionAsBdd.renameVariables(blockVariableSet, blockPrimeVariableSet), model.getColumnVariables()); // Pick a representative from each block. - partitionAsBdd = partitionAsBdd.existsAbstractRepresentative(model.getColumnVariables()); - partitionAsAdd = partitionAsBdd.template toAdd(); + auto representatives = InternalRepresentativeComputer(partition, model.getRowVariables(), model.getColumnVariables()).getRepresentatives(); + partitionAsBdd = representatives && partitionAsBdd.renameVariables(model.getColumnVariables(), model.getRowVariables()); + storm::dd::Add partitionAsAdd = partitionAsBdd.template toAdd(); - quotientTransitionMatrix = quotientTransitionMatrix.multiplyMatrix(partitionAsAdd.renameVariables(model.getColumnVariables(), model.getRowVariables()), model.getRowVariables()); + quotientTransitionMatrix = quotientTransitionMatrix.multiplyMatrix(partitionAsAdd, model.getRowVariables()); end = std::chrono::high_resolution_clock::now(); // Check quotient matrix for sanity. @@ -776,61 +961,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Cannot extract quotient for this model type."); } } - - template - std::shared_ptr> QuotientExtractor::extractQuotientUsingOriginalVariables(storm::models::symbolic::Model const& model, Partition const& partition, PreservationInformation const& preservationInformation) { - auto modelType = model.getType(); - - if (modelType == storm::models::ModelType::Dtmc || modelType == storm::models::ModelType::Ctmc) { - std::set blockVariableSet = {partition.getBlockVariable()}; - std::set blockPrimeVariableSet = {partition.getPrimedBlockVariable()}; - std::vector> blockMetaVariablePairs = {std::make_pair(partition.getBlockVariable(), partition.getPrimedBlockVariable())}; - - storm::dd::Add partitionAsAdd = partition.storedAsBdd() ? partition.asBdd().template toAdd() : partition.asAdd(); - storm::dd::Add quotientTransitionMatrix = model.getTransitionMatrix().multiplyMatrix(partitionAsAdd, model.getColumnVariables()); - quotientTransitionMatrix = quotientTransitionMatrix.renameVariables(blockVariableSet, model.getColumnVariables()); - partitionAsAdd = partitionAsAdd / partitionAsAdd.sumAbstract(model.getColumnVariables()); - quotientTransitionMatrix = quotientTransitionMatrix.multiplyMatrix(partitionAsAdd, model.getRowVariables()); - quotientTransitionMatrix = quotientTransitionMatrix.renameVariables(blockVariableSet, model.getRowVariables()); - storm::dd::Bdd quotientTransitionMatrixBdd = quotientTransitionMatrix.notZero(); - - // Check quotient matrix for sanity. - STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one()).isZero(), "Illegal entries in quotient matrix."); - STORM_LOG_ASSERT(quotientTransitionMatrix.sumAbstract(blockPrimeVariableSet).equalModuloPrecision(quotientTransitionMatrix.notZero().existsAbstract(blockPrimeVariableSet).template toAdd(), ValueType(1e-6)), "Illegal non-probabilistic matrix."); - - storm::dd::Bdd partitionAsBdd = partition.storedAsBdd() ? partition.asBdd() : partition.asAdd().notZero(); - storm::dd::Bdd partitionAsBddOverRowVariables = partitionAsBdd.renameVariables(model.getColumnVariables(), model.getRowVariables()); - storm::dd::Bdd reachableStates = partitionAsBdd.existsAbstract(model.getColumnVariables()).renameVariables(blockVariableSet, model.getRowVariables()); - storm::dd::Bdd initialStates = (model.getInitialStates() && partitionAsBdd.renameVariables(model.getColumnVariables(), model.getRowVariables())).existsAbstract(model.getRowVariables()).renameVariables(blockVariableSet, model.getRowVariables()); - storm::dd::Bdd deadlockStates = !quotientTransitionMatrixBdd.existsAbstract(model.getColumnVariables()) && reachableStates; - - std::map> preservedLabelBdds; - for (auto const& label : preservationInformation.getLabels()) { - preservedLabelBdds.emplace(label, (model.getStates(label) && partitionAsBddOverRowVariables).existsAbstract(model.getRowVariables())); - } - for (auto const& expression : preservationInformation.getExpressions()) { - std::stringstream stream; - stream << expression; - std::string expressionAsString = stream.str(); - auto it = preservedLabelBdds.find(expressionAsString); - if (it != preservedLabelBdds.end()) { - STORM_LOG_WARN("Duplicate label '" << expressionAsString << "', dropping second label definition."); - } else { - preservedLabelBdds.emplace(stream.str(), (model.getStates(expression) && partitionAsBddOverRowVariables).existsAbstract(model.getRowVariables())); - } - } - - if (modelType == storm::models::ModelType::Dtmc) { - return std::shared_ptr>(new storm::models::symbolic::Dtmc(model.getManager().asSharedPointer(), reachableStates, initialStates, deadlockStates, quotientTransitionMatrix, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs(), preservedLabelBdds, {})); - } else { - return std::shared_ptr>(new storm::models::symbolic::Ctmc(model.getManager().asSharedPointer(), reachableStates, initialStates, deadlockStates, quotientTransitionMatrix, blockVariableSet, blockPrimeVariableSet, blockMetaVariablePairs, preservedLabelBdds, {})); - } - } else { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Cannot exctract quotient for this model type."); - } - } - template class QuotientExtractor; template class QuotientExtractor; diff --git a/src/storm/storage/dd/bisimulation/QuotientExtractor.h b/src/storm/storage/dd/bisimulation/QuotientExtractor.h index dda06a3d3..dcce31659 100644 --- a/src/storm/storage/dd/bisimulation/QuotientExtractor.h +++ b/src/storm/storage/dd/bisimulation/QuotientExtractor.h @@ -28,7 +28,6 @@ namespace storm { std::shared_ptr> extractDdQuotient(storm::models::symbolic::Model const& model, Partition const& partition, PreservationInformation const& preservationInformation); std::shared_ptr> extractQuotientUsingBlockVariables(storm::models::symbolic::Model const& model, Partition const& partition, PreservationInformation const& preservationInformation); - std::shared_ptr> extractQuotientUsingOriginalVariables(storm::models::symbolic::Model const& model, Partition const& partition, PreservationInformation const& preservationInformation); bool useRepresentatives; storm::settings::modules::BisimulationSettings::QuotientFormat quotientFormat;