Merge branch 'monolithic-dft' of https://sselab.de/lab9/private/git/storm into monolithic-dft

Former-commit-id: c668dfce66

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;
+

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;
+

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;

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;
+

src/builder/ExplicitDFTModelBuilder.cpp

@@ -15,42 +15,10 @@ namespace storm {
         }
         
         template <typename ValueType>
-        ExplicitDFTModelBuilder<ValueType>::ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft) : mDft(dft), mStates(((mDft.stateVectorSize() / 64) + 1) * 64, INITIAL_BUCKETSIZE) {
+        ExplicitDFTModelBuilder<ValueType>::ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> 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<std::vector<size_t>> symmetriesTmp;
-            std::vector<size_t> vecB;
-            std::vector<size_t> vecC;
-            std::vector<size_t> 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<storm::storage::DFTStateGenerationInfo>(mDft.buildStateGenerationInfo(vecB, symmetriesTmp));
+            mStateGenerationInfo = std::make_shared<storm::storage::DFTStateGenerationInfo>(mDft.buildStateGenerationInfo(symmetries));
         }
 
 
@@ -166,7 +134,7 @@ namespace storm {
         bool ExplicitDFTModelBuilder<ValueType>::exploreStates(std::queue<DFTStatePointer>& stateQueue, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<uint_fast64_t>& markovianStates, std::vector<ValueType>& exitRates) {
 
             // TODO Matthias: set Markovian states directly as bitvector?
-            std::map<size_t, ValueType> outgoingTransitions;
+            std::map<uint_fast64_t, ValueType> 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<storm::storage::DFTState<ValueType>>(*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;
         }

src/builder/ExplicitDFTModelBuilder.h

@@ -1,13 +1,13 @@
 #ifndef EXPLICITDFTMODELBUILDER_H
 #define	EXPLICITDFTMODELBUILDER_H
 
-#include "../storage/dft/DFT.h"
-
 #include <src/models/sparse/StateLabeling.h>
 #include <src/models/sparse/StandardRewardModel.h>
 #include <src/models/sparse/Model.h>
 #include <src/storage/SparseMatrix.h>
 #include <src/storage/BitVectorHashMap.h>
+#include <src/storage/dft/DFT.h>
+#include <src/storage/dft/SymmetricUnits.h>
 #include <boost/container/flat_set.hpp>
 #include <boost/optional/optional.hpp>
 #include <stack>
@@ -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<ValueType> const& mDft;
             std::shared_ptr<storm::storage::DFTStateGenerationInfo> mStateGenerationInfo;
-            storm::storage::BitVectorHashMap<size_t> mStates;
+            storm::storage::BitVectorHashMap<uint_fast64_t> mStates;
+            std::vector<uint_fast64_t> mPseudoStatesMapping;
             size_t newIndex = 0;
             bool mergeFailedStates = true;
             size_t failedIndex = 0;
@@ -65,7 +65,7 @@ namespace storm {
                 std::set<std::string> beLabels = {};
             };
             
-            ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft);
+            ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTIndependentSymmetries const& symmetries);
 
             std::shared_ptr<storm::models::sparse::Model<ValueType>> buildModel(LabelOptions const& labelOpts);
 

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<ValueType> 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<ValueType> result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates, stateRewardValues, checkTask.isOnlyInitialStatesRelevantSet());
+            std::vector<ValueType> result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates,
+                                                                       [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) {
+                                                                           return rewardModel.getTotalRewardVector(numberOfRows, transitionMatrix, maybeStates);
+                                                                       },
+                                                                       checkTask.isOnlyInitialStatesRelevantSet());
 
             // Construct check result.
             std::unique_ptr<CheckResult> checkResult(new ExplicitQuantitativeCheckResult<ValueType>(result));
@@ -597,6 +599,17 @@ namespace storm {
 
         template<typename SparseDtmcModelType>
         std::vector<typename SparseDtmcModelType::ValueType> SparseDtmcEliminationModelChecker<SparseDtmcModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::vector<ValueType>& stateRewardValues, bool computeForInitialStatesOnly) {
+            return computeReachabilityRewards(probabilityMatrix, backwardTransitions, initialStates, targetStates,
+                                              [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) {
+                                                  std::vector<ValueType> result(numberOfRows);
+                                                  storm::utility::vector::selectVectorValues(result, maybeStates, stateRewardValues);
+                                                  return result;
+                                              },
+                                              computeForInitialStatesOnly);
+        }
+
+        template<typename SparseDtmcModelType>
+        std::vector<typename SparseDtmcModelType::ValueType> SparseDtmcEliminationModelChecker<SparseDtmcModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly) {
 
             uint_fast64_t numberOfStates = probabilityMatrix.getRowCount();
             
@@ -639,9 +652,9 @@ namespace storm {
                 storm::storage::SparseMatrix<ValueType> submatrixTransposed = submatrix.transpose();
                 
                 // Project the state reward vector to all maybe-states.
-                storm::storage::BitVector phiStates(numberOfStates, true);
+                std::vector<ValueType> stateRewardValues = totalStateRewardVectorGetter(submatrix.getRowCount(), probabilityMatrix, maybeStates);
 
-                std::vector<ValueType> subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, phiStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates));
+                std::vector<ValueType> subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, trueStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates));
                 storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, subresult);
             }
             

src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h

@@ -91,6 +91,8 @@ namespace storm {
             
             static std::vector<ValueType> computeLongRunValues(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& maybeStates, bool computeResultsForInitialStatesOnly, std::vector<ValueType>& stateValues);
             
+            static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly);
+
             static std::vector<ValueType> computeReachabilityValues(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType>& values, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, bool computeResultsForInitialStatesOnly, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& oneStepProbabilitiesToTarget);
             
             static std::shared_ptr<StatePriorityQueue<ValueType>> createStatePriorityQueue(boost::optional<std::vector<uint_fast64_t>> const& stateDistances, storm::storage::FlexibleSparseMatrix<ValueType> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& states);

src/storage/BitVectorHashMap.cpp

@@ -140,6 +140,11 @@ namespace storm {
             return findOrAddAndGetBucket(key, value).first;
         }
         
+        template<class ValueType, class Hash1, class Hash2>
+        void BitVectorHashMap<ValueType, Hash1, Hash2>::setOrAdd(storm::storage::BitVector const& key, ValueType const& value) {
+            setOrAddAndGetBucket(key, value);
+        }
+        
         template<class ValueType, class Hash1, class Hash2>
         std::pair<ValueType, std::size_t> BitVectorHashMap<ValueType, Hash1, Hash2>::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<class ValueType, class Hash1, class Hash2>
+        std::size_t BitVectorHashMap<ValueType, Hash1, Hash2>::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<bool, std::size_t, bool> flagBucketTuple = this->findBucketToInsert<true>(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<class ValueType, class Hash1, class Hash2>
         ValueType BitVectorHashMap<ValueType, Hash1, Hash2>::getValue(storm::storage::BitVector const& key) const {
             std::pair<bool, std::size_t> flagBucketPair = this->findBucket(key);

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<ValueType, std::size_t> 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.
              *

src/storage/SparseMatrix.cpp

@@ -220,6 +220,73 @@ namespace storm {
             return lastColumn;
         }
         
+        // Debug method for printing the current matrix
+        template<typename ValueType>
+        void print(std::vector<typename SparseMatrix<ValueType>::index_type> const& rowGroupIndices, std::vector<MatrixEntry<typename SparseMatrix<ValueType>::index_type, typename SparseMatrix<ValueType>::value_type>> const& columnsAndValues, std::vector<typename SparseMatrix<ValueType>::index_type> const& rowIndications) {
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::index_type endRows;
+            // Iterate over all row groups.
+            for (typename SparseMatrix<ValueType>::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<ValueType>::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<ValueType>::index_type pos = rowIndications[i]; pos < endRows; ++pos) {
+                        std::cout << "(" << columnsAndValues[pos].getColumn() << ": " << columnsAndValues[pos].getValue() << ") ";
+                    }
+                    std::cout << std::endl;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        bool SparseMatrixBuilder<ValueType>::replaceColumns(std::vector<index_type> 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<typename ValueType>
+        void SparseMatrixBuilder<ValueType>::fixColumns() {
+            // Sort columns per row
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::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<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                  return a.getColumn() < b.getColumn();
+                              });
+                    // Assert no equal elements
+                    assert(std::is_sorted(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                                          [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                              return a.getColumn() <= b.getColumn();
+                                          }));
+                }
+            }
+        }
+        
         template<typename ValueType>
         SparseMatrix<ValueType>::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) {
             // Intentionally left empty.

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<index_type> 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();
         };
         
         /*!

src/storage/dft/DFT.cpp

@@ -68,20 +68,77 @@ namespace storm {
         }
 
         template<typename ValueType>
-        DFTStateGenerationInfo DFT<ValueType>::buildStateGenerationInfo(std::vector<size_t> const& subTreeRoots, std::vector<std::vector<size_t>> const& symmetries) const {
+        DFTStateGenerationInfo DFT<ValueType>::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<size_t> visitQueue;
-            std::set<size_t> 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<size_t> 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<size_t> 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<DFTDependency<ValueType> 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<typename ValueType>
-        size_t DFT<ValueType>::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, std::set<size_t>& visited, size_t stateIndex) const {
+        size_t DFT<ValueType>::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<typename ValueType>
+        size_t DFT<ValueType>::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& 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<DFTGate<ValueType>>(element)->children()) {
+                if (mElements[id]->isGate()) {
+                    for (auto const& child : std::static_pointer_cast<DFTGate<ValueType>>(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 << "]";
                     }

src/storage/dft/DFT.h

@@ -42,6 +42,8 @@ namespace storm {
             std::map<size_t, size_t> mSpareActivationIndex; // id spare representative -> index in state
             std::vector<size_t> mIdToStateIndex; // id -> index first bit in state
             std::map<size_t, size_t> mRestrictedItems;
+            std::vector<std::pair<size_t, std::vector<size_t>>> 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<size_t>& 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<size_t> 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<size_t> const& subTreeRoots, std::vector<std::vector<size_t>> 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<size_t>& visitQueue, std::set<size_t>& visited, size_t stateIndex) const;
+            size_t performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, storm::storage::BitVector& visited, size_t stateIndex) const;
             
             size_t stateVectorSize() const {
                 return mStateVectorSize;

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<ValueType>& 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());

src/storage/dft/DFTState.cpp

@@ -242,6 +242,37 @@ namespace storm {
             }
             return false;
         }
+        
+        template<typename ValueType>
+        bool DFTState<ValueType>::orderBySymmetry() {
+            bool changed = false;
+            for (size_t pos = 0; pos < mStateGenerationInfo.getSymmetrySize(); ++pos) {
+                // Check each symmetry
+                size_t length = mStateGenerationInfo.getSymmetryLength(pos);
+                std::vector<size_t> 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<double>;

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) {

src/storm-dyftee.cpp

@@ -28,17 +28,19 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false)
     std::vector<std::shared_ptr<storm::logic::Formula>> parsedFormulas = storm::parseFormulasForExplicit(property);
     std::vector<std::shared_ptr<const storm::logic::Formula>> formulas(parsedFormulas.begin(), parsedFormulas.end());
     assert(formulas.size() == 1);
-    std::cout << "parsed formula." << std::endl;
+    std::cout << "Parsed formula." << std::endl;
+    std::map<size_t, std::vector<std::vector<size_t>>> 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<ValueType> builder(dft);
+    storm::builder::ExplicitDFTModelBuilder<ValueType> builder(dft, symmetries);
     typename storm::builder::ExplicitDFTModelBuilder<ValueType>::LabelOptions labeloptions; // TODO initialize this with the formula
     std::shared_ptr<storm::models::sparse::Model<ValueType>> 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<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas, storm::storage::BisimulationType::Weak)->template as<storm::models::sparse::Ctmc<ValueType>>();
     }
     
-    
     model->printModelInformationToStream(std::cout);
     
     // Model checking