Removed Duplicated code

src/storm/solver/TopologicalLinearEquationSolver.cpp

@@ -63,15 +63,6 @@ namespace storm {
                 createSortedSccDecomposition(needAdaptPrecision);
             }
             
-            //std::cout << "Sorted SCC decomposition: " << std::endl;
-            //for (auto const& scc : *this->sortedSccDecomposition) {
-                //std::cout << "SCC: ";
-              //  for (auto const& row : scc) {
-                    //std::cout << row << "  ";
-               // }
-                //std::cout << std::endl;
-            //}
-            
             // We do not need to adapt the precision if all SCCs are trivial (i.e., the system is acyclic)
             needAdaptPrecision = needAdaptPrecision && (this->sortedSccDecomposition->size() != this->getMatrixRowCount());
             
@@ -119,90 +110,12 @@ namespace storm {
         template<typename ValueType>
         void TopologicalLinearEquationSolver<ValueType>::createSortedSccDecomposition(bool needLongestChainSize) const {
             // Obtain the scc decomposition
-            auto sccDecomposition = storm::storage::StronglyConnectedComponentDecomposition<ValueType>(*this->A);
-            
-            // Get a mapping from matrix row to the corresponding scc
-            STORM_LOG_THROW(sccDecomposition.size() < std::numeric_limits<uint32_t>::max(), storm::exceptions::UnexpectedException, "The number of SCCs is too large.");
-            std::vector<uint32_t> sccIndices(this->A->getRowCount(), std::numeric_limits<uint32_t>::max());
-            uint32_t sccIndex = 0;
-            for (auto const& scc : sccDecomposition) {
-                for (auto const& row : scc) {
-                    sccIndices[row] = sccIndex;
-                }
-                ++sccIndex;
-            }
-            
-            // Prepare the resulting set of sorted sccs
-            this->sortedSccDecomposition = std::make_unique<std::vector<storm::storage::StronglyConnectedComponent>>();
-            std::vector<storm::storage::StronglyConnectedComponent>& sortedSCCs = *this->sortedSccDecomposition;
-            sortedSCCs.reserve(sccDecomposition.size());
-            
-            // Find a topological sort via DFS.
-            storm::storage::BitVector unsortedSCCs(sccDecomposition.size(), true);
-            std::vector<uint32_t> sccStack, chainSizes;
+            this->sortedSccDecomposition = std::make_unique<storm::storage::StronglyConnectedComponentDecomposition<ValueType>>(*this->A);
             if (needLongestChainSize) {
-                chainSizes.resize(sccDecomposition.size(), 1u);
-            }
-            uint32_t longestChainSize = 0;
-            uint32_t const token = std::numeric_limits<uint32_t>::max();
-            std::set<uint64_t> successorSCCs;
-
-            for (uint32_t firstUnsortedScc = 0; firstUnsortedScc < unsortedSCCs.size(); firstUnsortedScc = unsortedSCCs.getNextSetIndex(firstUnsortedScc + 1)) {
-                
-                sccStack.push_back(firstUnsortedScc);
-                while (!sccStack.empty()) {
-                    uint32_t currentSccIndex = sccStack.back();
-                    if (currentSccIndex != token) {
-                        // Check whether the SCC is still unprocessed
-                        if (unsortedSCCs.get(currentSccIndex)) {
-                            // Explore the successors of the scc.
-                            storm::storage::StronglyConnectedComponent const& currentScc = sccDecomposition.getBlock(currentSccIndex);
-                            // We first push a token on the stack in order to recognize later when all successors of this SCC have been explored already.
-                            sccStack.push_back(token);
-                            // Now add all successors that are not already sorted.
-                            // Successors should only be added once, so we first prepare a set of them and add them afterwards.
-                            successorSCCs.clear();
-                            for (auto const& row : currentScc) {
-                                for (auto const& entry : this->A->getRow(row)) {
-                                    auto const& successorSCC = sccIndices[entry.getColumn()];
-                                    if (successorSCC != currentSccIndex && unsortedSCCs.get(successorSCC)) {
-                                        successorSCCs.insert(successorSCC);
-                                    }
-                                }
-                            }
-                            sccStack.insert(sccStack.end(), successorSCCs.begin(), successorSCCs.end());
-                            
-                        }
-                    } else {
-                        // all successors of the current scc have already been explored.
-                        sccStack.pop_back(); // pop the token
-                        
-                        currentSccIndex = sccStack.back();
-                        storm::storage::StronglyConnectedComponent& scc = sccDecomposition.getBlock(currentSccIndex);
-                        
-                        // Compute the longest chain size for this scc
-                        if (needLongestChainSize) {
-                            uint32_t& currentChainSize = chainSizes[currentSccIndex];
-                            for (auto const& row : scc) {
-                                for (auto const& entry : this->A->getRow(row)) {
-                                    auto const& successorSCC = sccIndices[entry.getColumn()];
-                                    if (successorSCC != currentSccIndex) {
-                                        currentChainSize = std::max(currentChainSize, chainSizes[successorSCC] + 1);
-                                    }
-                                }
-                            }
-                            longestChainSize = std::max(longestChainSize, currentChainSize);
-                        }
-                        
-                        unsortedSCCs.set(currentSccIndex, false);
-                        sccStack.pop_back(); // pop the current scc index
-                        sortedSCCs.push_back(std::move(scc));
-                    }
-                }
-            }
-            
-            if (longestChainSize > 0) {
-                this->longestSccChainSize = longestChainSize;
+                this->longestSccChainSize = 0;
+                this->sortedSccDecomposition->sortTopologically(*this->A, &(this->longestSccChainSize.get()));
+            } else {
+                this->sortedSccDecomposition->sortTopologically(*this->A);
             }
         }
         

src/storm/solver/TopologicalLinearEquationSolver.h

@@ -67,7 +67,7 @@ namespace storm {
             NativeMultiplier<ValueType> multiplier;
 
             // cached auxiliary data
-            mutable std::unique_ptr<std::vector<storm::storage::StronglyConnectedComponent>> sortedSccDecomposition;
+            mutable std::unique_ptr<storm::storage::StronglyConnectedComponentDecomposition<ValueType>> sortedSccDecomposition;
             mutable boost::optional<uint64_t> longestSccChainSize;
             mutable std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> sccSolver;
         };

src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp

@@ -73,15 +73,6 @@ namespace storm {
                 createSortedSccDecomposition(needAdaptPrecision);
             }
             
-            //std::cout << "Sorted SCC decomposition: " << std::endl;
-            //for (auto const& scc : *this->sortedSccDecomposition) {
-                //std::cout << "SCC: ";
-              //  for (auto const& row : scc) {
-                    //std::cout << row << "  ";
-               // }
-                //std::cout << std::endl;
-            //}
-            
             // We do not need to adapt the precision if all SCCs are trivial (i.e., the system is acyclic)
             needAdaptPrecision = needAdaptPrecision && (this->sortedSccDecomposition->size() != this->A->getRowGroupCount());
             
@@ -149,90 +140,12 @@ namespace storm {
         template<typename ValueType>
         void TopologicalMinMaxLinearEquationSolver<ValueType>::createSortedSccDecomposition(bool needLongestChainSize) const {
             // Obtain the scc decomposition
-            auto sccDecomposition = storm::storage::StronglyConnectedComponentDecomposition<ValueType>(*this->A);
-            
-            // Get a mapping from matrix row to the corresponding scc
-            STORM_LOG_THROW(sccDecomposition.size() < std::numeric_limits<uint32_t>::max(), storm::exceptions::UnexpectedException, "The number of SCCs is too large.");
-            std::vector<uint32_t> sccIndices(this->A->getRowCount(), std::numeric_limits<uint32_t>::max());
-            uint32_t sccIndex = 0;
-            for (auto const& scc : sccDecomposition) {
-                for (auto const& row : scc) {
-                    sccIndices[row] = sccIndex;
-                }
-                ++sccIndex;
-            }
-            
-            // Prepare the resulting set of sorted sccs
-            this->sortedSccDecomposition = std::make_unique<std::vector<storm::storage::StronglyConnectedComponent>>();
-            std::vector<storm::storage::StronglyConnectedComponent>& sortedSCCs = *this->sortedSccDecomposition;
-            sortedSCCs.reserve(sccDecomposition.size());
-            
-            // Find a topological sort via DFS.
-            storm::storage::BitVector unsortedSCCs(sccDecomposition.size(), true);
-            std::vector<uint32_t> sccStack, chainSizes;
+            this->sortedSccDecomposition = std::make_unique<storm::storage::StronglyConnectedComponentDecomposition<ValueType>>(*this->A);
             if (needLongestChainSize) {
-                chainSizes.resize(sccDecomposition.size(), 1u);
-            }
-            uint32_t longestChainSize = 0;
-            uint32_t const token = std::numeric_limits<uint32_t>::max();
-            std::set<uint64_t> successorSCCs;
-
-            for (uint32_t firstUnsortedScc = 0; firstUnsortedScc < unsortedSCCs.size(); firstUnsortedScc = unsortedSCCs.getNextSetIndex(firstUnsortedScc + 1)) {
-                
-                sccStack.push_back(firstUnsortedScc);
-                while (!sccStack.empty()) {
-                    uint32_t currentSccIndex = sccStack.back();
-                    if (currentSccIndex != token) {
-                        // Check whether the SCC is still unprocessed
-                        if (unsortedSCCs.get(currentSccIndex)) {
-                            // Explore the successors of the scc.
-                            storm::storage::StronglyConnectedComponent const& currentScc = sccDecomposition.getBlock(currentSccIndex);
-                            // We first push a token on the stack in order to recognize later when all successors of this SCC have been explored already.
-                            sccStack.push_back(token);
-                            // Now add all successors that are not already sorted.
-                            // Successors should only be added once, so we first prepare a set of them and add them afterwards.
-                            successorSCCs.clear();
-                            for (auto const& group : currentScc) {
-                                for (auto const& entry : this->A->getRowGroup(group)) {
-                                    auto const& successorSCC = sccIndices[entry.getColumn()];
-                                    if (successorSCC != currentSccIndex && unsortedSCCs.get(successorSCC)) {
-                                        successorSCCs.insert(successorSCC);
-                                    }
-                                }
-                            }
-                            sccStack.insert(sccStack.end(), successorSCCs.begin(), successorSCCs.end());
-                            
-                        }
-                    } else {
-                        // all successors of the current scc have already been explored.
-                        sccStack.pop_back(); // pop the token
-                        
-                        currentSccIndex = sccStack.back();
-                        storm::storage::StronglyConnectedComponent& scc = sccDecomposition.getBlock(currentSccIndex);
-                        
-                        // Compute the longest chain size for this scc
-                        if (needLongestChainSize) {
-                            uint32_t& currentChainSize = chainSizes[currentSccIndex];
-                            for (auto const& row : scc) {
-                                for (auto const& entry : this->A->getRow(row)) {
-                                    auto const& successorSCC = sccIndices[entry.getColumn()];
-                                    if (successorSCC != currentSccIndex) {
-                                        currentChainSize = std::max(currentChainSize, chainSizes[successorSCC] + 1);
-                                    }
-                                }
-                            }
-                            longestChainSize = std::max(longestChainSize, currentChainSize);
-                        }
-                        
-                        unsortedSCCs.set(currentSccIndex, false);
-                        sccStack.pop_back(); // pop the current scc index
-                        sortedSCCs.push_back(std::move(scc));
-                    }
-                }
-            }
-            
-            if (longestChainSize > 0) {
-                this->longestSccChainSize = longestChainSize;
+                this->longestSccChainSize = 0;
+                this->sortedSccDecomposition->sortTopologically(*this->A, &(this->longestSccChainSize.get()));
+            } else {
+                this->sortedSccDecomposition->sortTopologically(*this->A);
             }
         }
         

src/storm/solver/TopologicalMinMaxLinearEquationSolver.h

@@ -21,7 +21,6 @@ namespace storm {
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
             
-
             virtual void clearCache() const override;
 
             virtual void repeatedMultiply(Environment const& env, OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1) const override;
@@ -54,7 +53,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> const* A;
             
             // cached auxiliary data
-            mutable std::unique_ptr<std::vector<storm::storage::StronglyConnectedComponent>> sortedSccDecomposition;
+            mutable std::unique_ptr<storm::storage::StronglyConnectedComponentDecomposition<ValueType>> sortedSccDecomposition;
             mutable boost::optional<uint64_t> longestSccChainSize;
             mutable std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> sccSolver;
             mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRowGroupVector; // A.rowGroupCount() entries

src/storm/storage/StronglyConnectedComponentDecomposition.cpp

@@ -2,6 +2,9 @@
 #include "storm/models/sparse/Model.h"
 #include "storm/models/sparse/StandardRewardModel.h"
 #include "storm/adapters/RationalFunctionAdapter.h"
+#include "storm/utility/macros.h"
+
+#include "storm/exceptions/UnexpectedException.h"
 
 namespace storm {
     namespace storage {
@@ -226,6 +229,91 @@ namespace storm {
             }
         }
         
+        template <typename ValueType>
+        void StronglyConnectedComponentDecomposition<ValueType>::sortTopologically(storm::storage::SparseMatrix<ValueType> const& transitions, uint64_t* longestChainSize) {
+            
+            // Get a mapping from state to the corresponding scc
+            STORM_LOG_THROW(this->size() < std::numeric_limits<uint32_t>::max(), storm::exceptions::UnexpectedException, "The number of SCCs is too large.");
+            std::vector<uint32_t> sccIndices(transitions.getRowGroupCount(), std::numeric_limits<uint32_t>::max());
+            uint32_t sccIndex = 0;
+            for (auto const& scc : *this) {
+                for (auto const& state : scc) {
+                    sccIndices[state] = sccIndex;
+                }
+                ++sccIndex;
+            }
+            
+            // Prepare the resulting set of sorted sccs
+            std::vector<storm::storage::StronglyConnectedComponent> sortedSCCs;
+            sortedSCCs.reserve(this->size());
+            
+            // Find a topological sort via DFS.
+            storm::storage::BitVector unsortedSCCs(this->size(), true);
+            std::vector<uint32_t> sccStack, chainSizes;
+            if (longestChainSize != nullptr) {
+                chainSizes.resize(this->size(), 1u);
+                *longestChainSize = 0;
+            }
+            uint32_t const token = std::numeric_limits<uint32_t>::max();
+            std::set<uint64_t> successorSCCs;
+
+            for (uint32_t firstUnsortedScc = 0; firstUnsortedScc < unsortedSCCs.size(); firstUnsortedScc = unsortedSCCs.getNextSetIndex(firstUnsortedScc + 1)) {
+                
+                sccStack.push_back(firstUnsortedScc);
+                while (!sccStack.empty()) {
+                    uint32_t currentSccIndex = sccStack.back();
+                    if (currentSccIndex != token) {
+                        // Check whether the SCC is still unprocessed
+                        if (unsortedSCCs.get(currentSccIndex)) {
+                            // Explore the successors of the scc.
+                            storm::storage::StronglyConnectedComponent const& currentScc = this->getBlock(currentSccIndex);
+                            // We first push a token on the stack in order to recognize later when all successors of this SCC have been explored already.
+                            sccStack.push_back(token);
+                            // Now add all successors that are not already sorted.
+                            // Successors should only be added once, so we first prepare a set of them and add them afterwards.
+                            successorSCCs.clear();
+                            for (auto const& state : currentScc) {
+                                for (auto const& entry : transitions.getRowGroup(state)) {
+                                    auto const& successorSCC = sccIndices[entry.getColumn()];
+                                    if (successorSCC != currentSccIndex && unsortedSCCs.get(successorSCC)) {
+                                        successorSCCs.insert(successorSCC);
+                                    }
+                                }
+                            }
+                            sccStack.insert(sccStack.end(), successorSCCs.begin(), successorSCCs.end());
+                            
+                        }
+                    } else {
+                        // all successors of the current scc have already been explored.
+                        sccStack.pop_back(); // pop the token
+                        
+                        currentSccIndex = sccStack.back();
+                        storm::storage::StronglyConnectedComponent& scc = this->getBlock(currentSccIndex);
+                        
+                        // Compute the longest chain size for this scc
+                        if (longestChainSize != nullptr) {
+                            uint32_t& currentChainSize = chainSizes[currentSccIndex];
+                            for (auto const& state : scc) {
+                                for (auto const& entry : transitions.getRowGroup(state)) {
+                                    auto const& successorSCC = sccIndices[entry.getColumn()];
+                                    if (successorSCC != currentSccIndex) {
+                                        currentChainSize = std::max(currentChainSize, chainSizes[successorSCC] + 1);
+                                    }
+                                }
+                            }
+                            *longestChainSize = std::max<uint64_t>(*longestChainSize, currentChainSize);
+                        }
+                        
+                        unsortedSCCs.set(currentSccIndex, false);
+                        sccStack.pop_back(); // pop the current scc index
+                        sortedSCCs.push_back(std::move(scc));
+                    }
+                }
+            }
+            this->blocks = std::move(sortedSCCs);
+        }
+
+        
         // Explicitly instantiate the SCC decomposition.
         template class StronglyConnectedComponentDecomposition<double>;
         template StronglyConnectedComponentDecomposition<double>::StronglyConnectedComponentDecomposition(storm::models::sparse::Model<double> const& model, bool dropNaiveSccs, bool onlyBottomSccs);

src/storm/storage/StronglyConnectedComponentDecomposition.h

@@ -131,6 +131,14 @@ namespace storm {
              */
             StronglyConnectedComponentDecomposition& operator=(StronglyConnectedComponentDecomposition&& other);
             
+            
+            /*!
+             * Sorts the SCCs topologically: The ith block can only reach states in block j<i
+             * @param longestChainSize if not nullptr, this value is set to the length m of the longest
+             * chain of SCCs B_1 B_2 ... B_m such that B_i can reach B_(i-1).
+             */
+            void sortTopologically(storm::storage::SparseMatrix<ValueType> const& transitions, uint64_t* longestChainSize = nullptr);
+            
         private:
             /*!
              * Performs the SCC decomposition of the given model. As a side-effect this fills the vector of