diff --git a/examples/dft/pdep_symmetry.dft b/examples/dft/pdep_symmetry.dft new file mode 100644 index 000000000..43224d790 --- /dev/null +++ b/examples/dft/pdep_symmetry.dft @@ -0,0 +1,11 @@ +toplevel "A"; +"A" and "B" "B'"; +"B" and "C" "D" "PDEP"; +"B'" and "C'" "D'" "PDEP'"; +"PDEP" pdep=0.6 "B" "C"; +"PDEP'" pdep=0.6 "B'" "C'"; +"C" lambda=0.5 dorm=0; +"D" lambda=0.5 dorm=0; +"C'" lambda=0.5 dorm=0; +"D'" lambda=0.5 dorm=0; + diff --git a/examples/dft/symmetry3.dft b/examples/dft/symmetry3.dft new file mode 100644 index 000000000..edda963f9 --- /dev/null +++ b/examples/dft/symmetry3.dft @@ -0,0 +1,12 @@ +toplevel "A"; +"A" and "B" "B'" "B''"; +"B" and "C" "D"; +"B'" and "C'" "D'"; +"B''" and "C''" "D''"; +"C" lambda=0.5 dorm=0; +"D" lambda=0.5 dorm=0; +"C'" lambda=0.5 dorm=0; +"D'" lambda=0.5 dorm=0; +"C''" lambda=0.5 dorm=0; +"D''" lambda=0.5 dorm=0; + diff --git a/examples/dft/symmetry4.dft b/examples/dft/symmetry4.dft new file mode 100644 index 000000000..f401fc241 --- /dev/null +++ b/examples/dft/symmetry4.dft @@ -0,0 +1,12 @@ +toplevel "A"; +"A" and "B" "B'" "C" "C'"; +"B" and "D" "E"; +"B'" and "D'" "E'"; +"C" or "F"; +"C'" or "F'"; +"D" lambda=0.5 dorm=0; +"E" lambda=0.5 dorm=0; +"D'" lambda=0.5 dorm=0; +"E'" lambda=0.5 dorm=0; +"F" lambda=0.5 dorm=0; +"F'" lambda=0.5 dorm=0; diff --git a/examples/dft/symmetry5.dft b/examples/dft/symmetry5.dft new file mode 100644 index 000000000..18f700bb2 --- /dev/null +++ b/examples/dft/symmetry5.dft @@ -0,0 +1,21 @@ +toplevel "A"; +"A" and "BA" "BB" "BC" "BD" "BE" "BF"; +"BA" and "CA" "DA"; +"BB" and "CB" "DB"; +"BC" and "CC" "DC"; +"BD" and "CD" "DD"; +"BE" and "CE" "DE"; +"BF" and "CF" "DF"; +"CA" lambda=0.5 dorm=0; +"DA" lambda=0.5 dorm=0; +"CB" lambda=0.5 dorm=0; +"DB" lambda=0.5 dorm=0; +"CC" lambda=0.5 dorm=0; +"DC" lambda=0.5 dorm=0; +"CD" lambda=0.5 dorm=0; +"DD" lambda=0.5 dorm=0; +"CE" lambda=0.5 dorm=0; +"DE" lambda=0.5 dorm=0; +"CF" lambda=0.5 dorm=0; +"DF" lambda=0.5 dorm=0; + diff --git a/src/builder/ExplicitDFTModelBuilder.cpp b/src/builder/ExplicitDFTModelBuilder.cpp index 1840ba40a..f14408190 100644 --- a/src/builder/ExplicitDFTModelBuilder.cpp +++ b/src/builder/ExplicitDFTModelBuilder.cpp @@ -15,42 +15,10 @@ namespace storm { } template - ExplicitDFTModelBuilder::ExplicitDFTModelBuilder(storm::storage::DFT const& dft) : mDft(dft), mStates(((mDft.stateVectorSize() / 64) + 1) * 64, INITIAL_BUCKETSIZE) { + ExplicitDFTModelBuilder::ExplicitDFTModelBuilder(storm::storage::DFT const& dft, storm::storage::DFTIndependentSymmetries const& symmetries) : mDft(dft), mStates(((mDft.stateVectorSize() / 64) + 1) * 64, INITIAL_BUCKETSIZE) { // stateVectorSize is bound for size of bitvector - // Find symmetries - // TODO activate - // Currently using hack to test - std::vector> symmetriesTmp; - std::vector vecB; - std::vector vecC; - std::vector vecD; - if (false) { - for (size_t i = 0; i < mDft.nrElements(); ++i) { - std::string name = mDft.getElement(i)->name(); - size_t id = mDft.getElement(i)->id(); - if (boost::starts_with(name, "B")) { - vecB.push_back(id); - } else if (boost::starts_with(name, "C")) { - vecC.push_back(id); - } else if (boost::starts_with(name, "D")) { - vecD.push_back(id); - } - } - symmetriesTmp.push_back(vecB); - symmetriesTmp.push_back(vecC); - symmetriesTmp.push_back(vecD); - std::cout << "Found the following symmetries:" << std::endl; - for (auto const& symmetry : symmetriesTmp) { - for (auto const& elem : symmetry) { - std::cout << elem << " -> "; - } - std::cout << std::endl; - } - } else { - vecB.push_back(mDft.getTopLevelIndex()); - } - mStateGenerationInfo = std::make_shared(mDft.buildStateGenerationInfo(vecB, symmetriesTmp)); + mStateGenerationInfo = std::make_shared(mDft.buildStateGenerationInfo(symmetries)); } @@ -166,7 +134,7 @@ namespace storm { bool ExplicitDFTModelBuilder::exploreStates(std::queue& stateQueue, storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, std::vector& markovianStates, std::vector& exitRates) { // TODO Matthias: set Markovian states directly as bitvector? - std::map outgoingTransitions; + std::map outgoingTransitions; size_t rowOffset = 0; // Captures number of non-deterministic choices bool deterministic = true; @@ -253,23 +221,51 @@ namespace storm { newState->updateFailableDependencies(nextBE->id()); bool dftFailed = newState->hasFailed(mDft.getTopLevelIndex()); - size_t newStateId; + uint_fast64_t newStateId; if(dftFailed && mergeFailedStates) { newStateId = failedIndex; } else { - if (mStates.contains(newState->status())) { + // Order state by symmetry + STORM_LOG_TRACE("Check for symmetry: " << mDft.getStateString(newState)); + bool changed = newState->orderBySymmetry(); + STORM_LOG_TRACE("State " << (changed ? "changed to " : "did not change") << (changed ? mDft.getStateString(newState) : "")); + + // Check if state already exists + if (mStates.contains(newState->status())) { // State already exists newStateId = mStates.getValue(newState->status()); - STORM_LOG_TRACE("State " << mDft.getStateString(newState) << " already exists"); + STORM_LOG_TRACE("State " << mDft.getStateString(newState) << " with id " << newStateId << " already exists"); + + // Check if possible pseudo state can be created now + if (!changed && newStateId >= OFFSET_PSEUDO_STATE) { + newStateId = newStateId - OFFSET_PSEUDO_STATE; + assert(newStateId < mPseudoStatesMapping.size()); + if (mPseudoStatesMapping[newStateId] == 0) { + // Create pseudo state now + newState->setId(newIndex++); + mPseudoStatesMapping[newStateId] = newState->getId(); + newStateId = newState->getId(); + mStates.setOrAdd(newState->status(), newStateId); + stateQueue.push(newState); + STORM_LOG_TRACE("Now create state " << mDft.getStateString(newState) << " with id " << newStateId); + } + } } else { // New state - newState->setId(newIndex++); - newStateId = mStates.findOrAdd(newState->status(), newState->getId()); - STORM_LOG_TRACE("New state " << mDft.getStateString(newState)); - - // Add state to search queue - stateQueue.push(newState); - } + if (changed) { + // Remember to create state later on + newState->setId(mPseudoStatesMapping.size() + OFFSET_PSEUDO_STATE); + mPseudoStatesMapping.push_back(0); + newStateId = mStates.findOrAdd(newState->status(), newState->getId()); + STORM_LOG_TRACE("New state for later creation: " << mDft.getStateString(newState)); + } else { + // Create new state + newState->setId(newIndex++); + newStateId = mStates.findOrAdd(newState->status(), newState->getId()); + STORM_LOG_TRACE("New state: " << mDft.getStateString(newState)); + stateQueue.push(newState); + } + } } // Set transitions @@ -280,10 +276,11 @@ namespace storm { STORM_LOG_TRACE("Added transition from " << state->getId() << " to " << newStateId << " with probability " << dependency->probability()); if (!storm::utility::isOne(dependency->probability())) { + // TODO Matthias: use symmetry as well // Add transition to state where dependency was unsuccessful DFTStatePointer unsuccessfulState = std::make_shared>(*state); unsuccessfulState->letDependencyBeUnsuccessful(smallest-1); - size_t unsuccessfulStateId; + uint_fast64_t unsuccessfulStateId; if (mStates.contains(unsuccessfulState->status())) { // Unsuccessful state already exists unsuccessfulStateId = mStates.getValue(unsuccessfulState->status()); @@ -335,6 +332,38 @@ namespace storm { if (hasDependencies) { rowOffset--; // One increment to many } else { + // Try to merge pseudo states with their instantiation + // TODO Matthias: improve? + auto it = outgoingTransitions.begin(); + while (it != outgoingTransitions.end()) { + if (it->first >= OFFSET_PSEUDO_STATE) { + uint_fast64_t newId = it->first - OFFSET_PSEUDO_STATE; + assert(newId < mPseudoStatesMapping.size()); + if (mPseudoStatesMapping[newId] > 0) { + // State exists already + newId = mPseudoStatesMapping[newId]; + auto itFind = outgoingTransitions.find(newId); + if (itFind != outgoingTransitions.end()) { + // Add probability from pseudo state to instantiation + itFind->second += it->second; + STORM_LOG_TRACE("Merged pseudo state " << newId << " adding rate " << it->second << " to total rate of " << itFind->second); + } else { + // Only change id + outgoingTransitions.emplace(newId, it->second); + STORM_LOG_TRACE("Instantiated pseudo state " << newId << " with rate " << it->second); + } + // Remove pseudo state + auto itErase = it; + ++it; + outgoingTransitions.erase(itErase); + } else { + ++it; + } + } else { + ++it; + } + } + // Add all probability transitions for (auto it = outgoingTransitions.begin(); it != outgoingTransitions.end(); ++it) { @@ -348,6 +377,11 @@ namespace storm { assert(exitRates.size()-1 == state->getId()); } // end while queue + assert(std::find(mPseudoStatesMapping.begin(), mPseudoStatesMapping.end(), 0) == mPseudoStatesMapping.end()); + + // Replace pseudo states in matrix + transitionMatrixBuilder.replaceColumns(mPseudoStatesMapping, OFFSET_PSEUDO_STATE); + assert(!deterministic || rowOffset == 0); return deterministic; } diff --git a/src/builder/ExplicitDFTModelBuilder.h b/src/builder/ExplicitDFTModelBuilder.h index 42ca2755d..b9008f91e 100644 --- a/src/builder/ExplicitDFTModelBuilder.h +++ b/src/builder/ExplicitDFTModelBuilder.h @@ -1,13 +1,13 @@ #ifndef EXPLICITDFTMODELBUILDER_H #define EXPLICITDFTMODELBUILDER_H -#include "../storage/dft/DFT.h" - #include #include #include #include #include +#include +#include #include #include #include @@ -47,12 +47,12 @@ namespace storm { }; const size_t INITIAL_BUCKETSIZE = 20000; - - + const uint_fast64_t OFFSET_PSEUDO_STATE = UINT_FAST64_MAX / 2; storm::storage::DFT const& mDft; std::shared_ptr mStateGenerationInfo; - storm::storage::BitVectorHashMap mStates; + storm::storage::BitVectorHashMap mStates; + std::vector mPseudoStatesMapping; size_t newIndex = 0; bool mergeFailedStates = true; size_t failedIndex = 0; @@ -65,7 +65,7 @@ namespace storm { std::set beLabels = {}; }; - ExplicitDFTModelBuilder(storm::storage::DFT const& dft); + ExplicitDFTModelBuilder(storm::storage::DFT const& dft, storm::storage::DFTIndependentSymmetries const& symmetries); std::shared_ptr> buildModel(LabelOptions const& labelOpts); diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index 2ea8fa459..5425044b2 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -584,11 +584,13 @@ namespace storm { // Do some sanity checks to establish some required properties. RewardModelType const& rewardModel = this->getModel().getRewardModel(checkTask.isRewardModelSet() ? checkTask.getRewardModel() : ""); - // TODO use maybeStates instead of all states - std::vector stateRewardValues = rewardModel.getTotalRewardVector(trueStates.getNumberOfSetBits(), this->getModel().getTransitionMatrix(), trueStates); STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::IllegalArgumentException, "Input model does not have a reward model."); - std::vector result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates, stateRewardValues, checkTask.isOnlyInitialStatesRelevantSet()); + std::vector result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates, + [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates) { + return rewardModel.getTotalRewardVector(numberOfRows, transitionMatrix, maybeStates); + }, + checkTask.isOnlyInitialStatesRelevantSet()); // Construct check result. std::unique_ptr checkResult(new ExplicitQuantitativeCheckResult(result)); @@ -597,6 +599,17 @@ namespace storm { template std::vector SparseDtmcEliminationModelChecker::computeReachabilityRewards(storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::vector& stateRewardValues, bool computeForInitialStatesOnly) { + return computeReachabilityRewards(probabilityMatrix, backwardTransitions, initialStates, targetStates, + [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates) { + std::vector result(numberOfRows); + storm::utility::vector::selectVectorValues(result, maybeStates, stateRewardValues); + return result; + }, + computeForInitialStatesOnly); + } + + template + std::vector SparseDtmcEliminationModelChecker::computeReachabilityRewards(storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly) { uint_fast64_t numberOfStates = probabilityMatrix.getRowCount(); @@ -639,9 +652,9 @@ namespace storm { storm::storage::SparseMatrix submatrixTransposed = submatrix.transpose(); // Project the state reward vector to all maybe-states. - storm::storage::BitVector phiStates(numberOfStates, true); + std::vector stateRewardValues = totalStateRewardVectorGetter(submatrix.getRowCount(), probabilityMatrix, maybeStates); - std::vector subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, phiStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates)); + std::vector subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, trueStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates)); storm::utility::vector::setVectorValues(result, maybeStates, subresult); } diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h index 76e939d66..33b61a2c4 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h @@ -91,6 +91,8 @@ namespace storm { static std::vector computeLongRunValues(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& maybeStates, bool computeResultsForInitialStatesOnly, std::vector& stateValues); + static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly); + static std::vector computeReachabilityValues(storm::storage::SparseMatrix const& transitionMatrix, std::vector& values, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, bool computeResultsForInitialStatesOnly, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& oneStepProbabilitiesToTarget); static std::shared_ptr> createStatePriorityQueue(boost::optional> const& stateDistances, storm::storage::FlexibleSparseMatrix const& transitionMatrix, storm::storage::FlexibleSparseMatrix const& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector const& states); diff --git a/src/storage/BitVectorHashMap.cpp b/src/storage/BitVectorHashMap.cpp index 4d6238350..08875af3b 100644 --- a/src/storage/BitVectorHashMap.cpp +++ b/src/storage/BitVectorHashMap.cpp @@ -140,6 +140,11 @@ namespace storm { return findOrAddAndGetBucket(key, value).first; } + template + void BitVectorHashMap::setOrAdd(storm::storage::BitVector const& key, ValueType const& value) { + setOrAddAndGetBucket(key, value); + } + template std::pair BitVectorHashMap::findOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value) { // If the load of the map is too high, we increase the size. @@ -161,6 +166,25 @@ namespace storm { } } + template + std::size_t BitVectorHashMap::setOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value) { + // If the load of the map is too high, we increase the size. + if (numberOfElements >= loadFactor * *currentSizeIterator) { + this->increaseSize(); + } + + std::tuple flagBucketTuple = this->findBucketToInsert(key); + STORM_LOG_ASSERT(!std::get<2>(flagBucketTuple), "Failed to find bucket for insertion."); + if (!std::get<0>(flagBucketTuple)) { + // Insert the new bits into the bucket. + buckets.set(std::get<1>(flagBucketTuple) * bucketSize, key); + occupied.set(std::get<1>(flagBucketTuple)); + ++numberOfElements; + } + values[std::get<1>(flagBucketTuple)] = value; + return std::get<1>(flagBucketTuple); + } + template ValueType BitVectorHashMap::getValue(storm::storage::BitVector const& key) const { std::pair flagBucketPair = this->findBucket(key); diff --git a/src/storage/BitVectorHashMap.h b/src/storage/BitVectorHashMap.h index d50dac59e..18ca52186 100644 --- a/src/storage/BitVectorHashMap.h +++ b/src/storage/BitVectorHashMap.h @@ -66,6 +66,15 @@ namespace storm { * @return The found value if the key is already contained in the map and the provided new value otherwise. */ ValueType findOrAdd(storm::storage::BitVector const& key, ValueType const& value); + + /*! + * Sets the given key value pain in the map. If the key is found in the map, the corresponding value is + * overwritten with the given value. Otherwise, the key is inserted with the given value. + * + * @param key The key to search or insert. + * @param value The value to set. + */ + void setOrAdd(storm::storage::BitVector const& key, ValueType const& value); /*! * Searches for the given key in the map. If it is found, the mapped-to value is returned. Otherwise, the @@ -79,6 +88,16 @@ namespace storm { */ std::pair findOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value); + /*! + * Sets the given key value pain in the map. If the key is found in the map, the corresponding value is + * overwritten with the given value. Otherwise, the key is inserted with the given value. + * + * @param key The key to search or insert. + * @param value The value to set. + * @return The index of the bucket into which the key was inserted. + */ + std::size_t setOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value); + /*! * Retrieves the key stored in the given bucket (if any) and the value it is mapped to. * diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index d41e2ef2c..fa95f3930 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -220,6 +220,73 @@ namespace storm { return lastColumn; } + // Debug method for printing the current matrix + template + void print(std::vector::index_type> const& rowGroupIndices, std::vector::index_type, typename SparseMatrix::value_type>> const& columnsAndValues, std::vector::index_type> const& rowIndications) { + typename SparseMatrix::index_type endGroups; + typename SparseMatrix::index_type endRows; + // Iterate over all row groups. + for (typename SparseMatrix::index_type group = 0; group < rowGroupIndices.size(); ++group) { + std::cout << "\t---- group " << group << "/" << (rowGroupIndices.size() - 1) << " ---- " << std::endl; + endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size(); + // Iterate over all rows in a row group + for (typename SparseMatrix::index_type i = rowGroupIndices[group]; i < endGroups; ++i) { + endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size(); + // Print the actual row. + std::cout << "Row " << i << " (" << rowIndications[i] << " - " << endRows << ")" << ": "; + for (typename SparseMatrix::index_type pos = rowIndications[i]; pos < endRows; ++pos) { + std::cout << "(" << columnsAndValues[pos].getColumn() << ": " << columnsAndValues[pos].getValue() << ") "; + } + std::cout << std::endl; + } + } + } + + template + bool SparseMatrixBuilder::replaceColumns(std::vector const& replacements, index_type offset) { + bool changed = false; + index_type maxColumn = 0; + for (auto& elem : columnsAndValues) { + if (elem.getColumn() >= offset) { + elem.setColumn(replacements[elem.getColumn() - offset]); + changed = true; + } + maxColumn = std::max(maxColumn, elem.getColumn()); + } + assert(changed || highestColumn == maxColumn); + highestColumn = maxColumn; + assert(changed || lastColumn == columnsAndValues[columnsAndValues.size() - 1].getColumn()); + lastColumn = columnsAndValues[columnsAndValues.size() - 1].getColumn(); + + if (changed) { + fixColumns(); + } + return changed; + } + + template + void SparseMatrixBuilder::fixColumns() { + // Sort columns per row + typename SparseMatrix::index_type endGroups; + typename SparseMatrix::index_type endRows; + for (index_type group = 0; group < rowGroupIndices.size(); ++group) { + endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size(); + for (index_type i = rowGroupIndices[group]; i < endGroups; ++i) { + endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size(); + // Sort the row + std::sort(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows, + [](MatrixEntry const& a, MatrixEntry const& b) { + return a.getColumn() < b.getColumn(); + }); + // Assert no equal elements + assert(std::is_sorted(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows, + [](MatrixEntry const& a, MatrixEntry const& b) { + return a.getColumn() <= b.getColumn(); + })); + } + } + } + template SparseMatrix::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) { // Intentionally left empty. diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index f9d9916c4..4892dbe32 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -217,6 +217,17 @@ namespace storm { */ index_type getLastColumn() const; + /*! + * Replaces all columns with id > offset according to replacements. + * Every state with id offset+i is replaced by the id in replacements[i]. + * Afterwards the columns are sorted. + * + * @param replacements Mapping indicating the replacements from offset+i -> value of i. + * @param offset Offset to add to each id in vector index. + * @return True if replacement took place, False if nothing changed. + */ + bool replaceColumns(std::vector const& replacements, index_type offset); + private: // A flag indicating whether a row count was set upon construction. bool initialRowCountSet; @@ -277,6 +288,11 @@ namespace storm { // Stores the currently active row group. This is used for correctly constructing the row grouping of the // matrix. index_type currentRowGroup; + + /*! + * Fixes the matrix by sorting the columns to gain increasing order again. + */ + void fixColumns(); }; /*! diff --git a/src/storage/dft/DFT.cpp b/src/storage/dft/DFT.cpp index 3fbcd5ee1..eef83928d 100644 --- a/src/storage/dft/DFT.cpp +++ b/src/storage/dft/DFT.cpp @@ -68,20 +68,77 @@ namespace storm { } template - DFTStateGenerationInfo DFT::buildStateGenerationInfo(std::vector const& subTreeRoots, std::vector> const& symmetries) const { + DFTStateGenerationInfo DFT::buildStateGenerationInfo(storm::storage::DFTIndependentSymmetries const& symmetries) const { // Use symmetry // Collect all elements in the first subtree // TODO make recursive to use for nested subtrees - + DFTStateGenerationInfo generationInfo(nrElements()); // Perform DFS and insert all elements of subtree sequentially size_t stateIndex = 0; std::queue visitQueue; - std::set visited; - visitQueue.push(subTreeRoots[0]); - stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex); - + storm::storage::BitVector visited(nrElements(), false); + + if (symmetries.groups.empty()) { + // Perform DFS for whole tree + visitQueue.push(mTopLevelIndex); + stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex); + } else { + for (auto const& symmetryGroup : symmetries.groups) { + assert(!symmetryGroup.second.empty()); + + // Perform DFS for first subtree of each symmetry + visitQueue.push(symmetryGroup.first); + size_t groupIndex = stateIndex; + stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex); + size_t offset = stateIndex - groupIndex; + + // Mirror symmetries + size_t noSymmetricElements = symmetryGroup.second.front().size(); + assert(noSymmetricElements > 1); + + for (std::vector symmetricElements : symmetryGroup.second) { + assert(symmetricElements.size() == noSymmetricElements); + + // Initialize for original element + size_t originalElement = symmetricElements[0]; + size_t index = generationInfo.getStateIndex(originalElement); + size_t activationIndex = isRepresentative(originalElement) ? generationInfo.getSpareActivationIndex(originalElement) : 0; + size_t usageIndex = mElements[originalElement]->isSpareGate() ? generationInfo.getSpareUsageIndex(originalElement) : 0; + + // Mirror symmetry for each element + for (size_t i = 1; i < symmetricElements.size(); ++i) { + size_t symmetricElement = symmetricElements[i]; + visited.set(symmetricElement); + + generationInfo.addStateIndex(symmetricElement, index + offset * i); + stateIndex += 2; + + assert((activationIndex > 0) == isRepresentative(symmetricElement)); + if (activationIndex > 0) { + generationInfo.addSpareActivationIndex(symmetricElement, activationIndex + offset * i); + ++stateIndex; + } + + assert((usageIndex > 0) == mElements[symmetricElement]->isSpareGate()); + if (usageIndex > 0) { + generationInfo.addSpareUsageIndex(symmetricElement, usageIndex + offset * i); + stateIndex += generationInfo.usageInfoBits(); + } + } + } + + // Store starting indices of symmetry groups + std::vector symmetryIndices; + for (size_t i = 0; i < noSymmetricElements; ++i) { + symmetryIndices.push_back(groupIndex + i * offset); + } + generationInfo.addSymmetry(offset, symmetryIndices); + } + } + + // TODO symmetries in dependencies? // Consider dependencies for (size_t idDependency : getDependencies()) { std::shared_ptr const> dependency = getDependency(idDependency); @@ -90,47 +147,61 @@ namespace storm { visitQueue.push(dependency->dependentEvent()->id()); } stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex); - - assert(stateIndex = mStateVectorSize); + + // Visit all remaining states + for (size_t i = 0; i < visited.size(); ++i) { + if (!visited[i]) { + visitQueue.push(i); + stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex); + } + } STORM_LOG_TRACE(generationInfo); + assert(stateIndex == mStateVectorSize); + assert(visited.full()); return generationInfo; } template - size_t DFT::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue& visitQueue, std::set& visited, size_t stateIndex) const { + size_t DFT::generateStateInfo(DFTStateGenerationInfo& generationInfo, size_t id, storm::storage::BitVector& visited, size_t stateIndex) const { + assert(!visited[id]); + visited.set(id); + + // Reserve bits for element + generationInfo.addStateIndex(id, stateIndex); + stateIndex += 2; + + if (isRepresentative(id)) { + generationInfo.addSpareActivationIndex(id, stateIndex); + ++stateIndex; + } + + if (mElements[id]->isSpareGate()) { + generationInfo.addSpareUsageIndex(id, stateIndex); + stateIndex += generationInfo.usageInfoBits(); + } + + return stateIndex; + } + + template + size_t DFT::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue& visitQueue, storm::storage::BitVector& visited, size_t stateIndex) const { while (!visitQueue.empty()) { size_t id = visitQueue.front(); visitQueue.pop(); - if (visited.count(id) == 1) { + if (visited[id]) { // Already visited continue; } - visited.insert(id); - DFTElementPointer element = mElements[id]; + stateIndex = generateStateInfo(generationInfo, id, visited, stateIndex); // Insert children - if (element->isGate()) { - for (auto const& child : std::static_pointer_cast>(element)->children()) { + if (mElements[id]->isGate()) { + for (auto const& child : std::static_pointer_cast>(mElements[id])->children()) { visitQueue.push(child->id()); } } - - // Reserve bits for element - generationInfo.addStateIndex(id, stateIndex); - stateIndex += 2; - - if (isRepresentative(id)) { - generationInfo.addSpareActivationIndex(id, stateIndex); - ++stateIndex; - } - - if (element->isSpareGate()) { - generationInfo.addSpareUsageIndex(id, stateIndex); - stateIndex += generationInfo.usageInfoBits(); - } - } return stateIndex; } @@ -192,10 +263,10 @@ namespace storm { stream << "\t** " << storm::storage::toChar(state->getElementState(elem->id())); if(elem->isSpareGate()) { size_t useId = state->uses(elem->id()); - if(state->isActive(useId)) { - stream << " actively "; + if(useId == elem->id() || state->isActive(useId)) { + stream << "actively "; } - stream << " using " << useId; + stream << "using " << useId; } } stream << std::endl; @@ -216,8 +287,8 @@ namespace storm { if(elem->isSpareGate()) { stream << "["; size_t useId = state->uses(elem->id()); - if(state->isActive(useId)) { - stream << " actively "; + if(useId == elem->id() || state->isActive(useId)) { + stream << "actively "; } stream << "using " << useId << "]"; } diff --git a/src/storage/dft/DFT.h b/src/storage/dft/DFT.h index f0e08b31a..e8ac01f19 100644 --- a/src/storage/dft/DFT.h +++ b/src/storage/dft/DFT.h @@ -42,6 +42,8 @@ namespace storm { std::map mSpareActivationIndex; // id spare representative -> index in state std::vector mIdToStateIndex; // id -> index first bit in state std::map mRestrictedItems; + std::vector>> mSymmetries; // pair (lenght of symmetry group, vector indicating the starting points of the symmetry groups) + public: DFTStateGenerationInfo(size_t nrElements) : mUsageInfoBits(nrElements > 1 ? storm::utility::math::uint64_log2(nrElements-1) + 1 : 1), mIdToStateIndex(nrElements) { @@ -79,6 +81,24 @@ namespace storm { return mSpareActivationIndex.at(id); } + size_t addSymmetry(size_t lenght, std::vector& startingIndices) { + mSymmetries.push_back(std::make_pair(lenght, startingIndices)); + } + + size_t getSymmetrySize() const { + return mSymmetries.size(); + } + + size_t getSymmetryLength(size_t pos) const { + assert(pos < mSymmetries.size()); + return mSymmetries[pos].first; + } + + std::vector const& getSymmetryIndices(size_t pos) const { + assert(pos < mSymmetries.size()); + return mSymmetries[pos].second; + } + friend std::ostream& operator<<(std::ostream& os, DFTStateGenerationInfo const& info) { os << "Id to state index:" << std::endl; for (size_t id = 0; id < info.mIdToStateIndex.size(); ++id) { @@ -92,6 +112,14 @@ namespace storm { for (auto pair : info.mSpareActivationIndex) { os << pair.first << " -> " << pair.second << std::endl; } + os << "Symmetries:" << std::endl; + for (auto pair : info.mSymmetries) { + os << "Length: " << pair.first << ", starting indices: "; + for (size_t index : pair.second) { + os << index << ", "; + } + os << std::endl; + } return os; } @@ -126,9 +154,11 @@ namespace storm { public: DFT(DFTElementVector const& elements, DFTElementPointer const& tle); - DFTStateGenerationInfo buildStateGenerationInfo(std::vector const& subTreeRoots, std::vector> const& symmetries) const; + DFTStateGenerationInfo buildStateGenerationInfo(storm::storage::DFTIndependentSymmetries const& symmetries) const; + + size_t generateStateInfo(DFTStateGenerationInfo& generationInfo, size_t id, storm::storage::BitVector& visited, size_t stateIndex) const; - size_t performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue& visitQueue, std::set& visited, size_t stateIndex) const; + size_t performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue& visitQueue, storm::storage::BitVector& visited, size_t stateIndex) const; size_t stateVectorSize() const { return mStateVectorSize; diff --git a/src/storage/dft/DFTElements.h b/src/storage/dft/DFTElements.h index 00688f36e..7867f399f 100644 --- a/src/storage/dft/DFTElements.h +++ b/src/storage/dft/DFTElements.h @@ -446,12 +446,12 @@ namespace storm { } /** - * Finish failed/failsafe spare gate by activating the children and setting the useIndex to zero. + * Finish failed/failsafe spare gate by activating the children and setting the useIndex to the spare id. * This prevents multiple fail states with different usages or activations. * @param state The current state. */ void finalizeSpare(DFTState& state) const { - state.setUses(this->mId, 0); + state.setUses(this->mId, this->mId); for (auto child : this->children()) { if (!state.isActive(child->id())) { state.activate(child->id()); diff --git a/src/storage/dft/DFTState.cpp b/src/storage/dft/DFTState.cpp index 7f34e03ef..be5c33060 100644 --- a/src/storage/dft/DFTState.cpp +++ b/src/storage/dft/DFTState.cpp @@ -242,6 +242,37 @@ namespace storm { } return false; } + + template + bool DFTState::orderBySymmetry() { + bool changed = false; + for (size_t pos = 0; pos < mStateGenerationInfo.getSymmetrySize(); ++pos) { + // Check each symmetry + size_t length = mStateGenerationInfo.getSymmetryLength(pos); + std::vector symmetryIndices = mStateGenerationInfo.getSymmetryIndices(pos); + // Sort symmetry group in decreasing order by bubble sort + // TODO use better algorithm? + size_t tmp, elem1, elem2; + size_t n = symmetryIndices.size(); + do { + tmp = 0; + for (size_t i = 1; i < n; ++i) { + elem1 = mStatus.getAsInt(symmetryIndices[i-1], length); + elem2 = mStatus.getAsInt(symmetryIndices[i], length); + if (elem1 < elem2) { + // Swap elements + mStatus.setFromInt(symmetryIndices[i-1], length, elem2); + mStatus.setFromInt(symmetryIndices[i], length, elem1); + tmp = i; + changed = true; + } + } + n = tmp; + } while (n > 0); + } + return changed; + } + // Explicitly instantiate the class. template class DFTState; diff --git a/src/storage/dft/DFTState.h b/src/storage/dft/DFTState.h index 7a26e10ab..eccf7f01d 100644 --- a/src/storage/dft/DFTState.h +++ b/src/storage/dft/DFTState.h @@ -164,6 +164,12 @@ namespace storm { */ void letDependencyBeUnsuccessful(size_t index = 0); + /** + * Order the state in decreasing order using the symmetries. + * @return True, if elements were swapped, false if nothing changed. + */ + bool orderBySymmetry(); + std::string getCurrentlyFailableString() const { std::stringstream stream; if (nrFailableDependencies() > 0) { diff --git a/src/storm-dyftee.cpp b/src/storm-dyftee.cpp index aa2e0879f..ba50d90f3 100644 --- a/src/storm-dyftee.cpp +++ b/src/storm-dyftee.cpp @@ -28,17 +28,19 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false) std::vector> parsedFormulas = storm::parseFormulasForExplicit(property); std::vector> formulas(parsedFormulas.begin(), parsedFormulas.end()); assert(formulas.size() == 1); - std::cout << "parsed formula." << std::endl; + std::cout << "Parsed formula." << std::endl; + std::map>> emptySymmetry; + storm::storage::DFTIndependentSymmetries symmetries(emptySymmetry); if(symred) { std::cout << dft.getElementsString() << std::endl; auto colouring = dft.colourDFT(); - storm::storage::DFTIndependentSymmetries symmetries = dft.findSymmetries(colouring); - std::cout << symmetries; + symmetries = dft.findSymmetries(colouring); + std::cout << "Symmetries: " << symmetries << std::endl; } // Building Markov Automaton std::cout << "Building Model..." << std::endl; - storm::builder::ExplicitDFTModelBuilder builder(dft); + storm::builder::ExplicitDFTModelBuilder builder(dft, symmetries); typename storm::builder::ExplicitDFTModelBuilder::LabelOptions labeloptions; // TODO initialize this with the formula std::shared_ptr> model = builder.buildModel(labeloptions); std::cout << "Built Model" << std::endl; @@ -52,7 +54,6 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false) model = storm::performDeterministicSparseBisimulationMinimization>(model->template as>(), formulas, storm::storage::BisimulationType::Weak)->template as>(); } - model->printModelInformationToStream(std::cout); // Model checking