REA fully implemented; needs testing

Former-commit-id: 9795a24835 [formerly fc732962dd]
Former-commit-id: 0ae2abacd1

src/test/utility/GraphTest.cpp

@@ -276,6 +276,12 @@ TEST(GraphTest, kshortest) {
 
     storm::utility::shortestPaths::ShortestPathsGenerator<double> shortestPathsGenerator(model);
 
+    storm::storage::sparse::state_type exampleState = 50;
+
+    for (int i = 1; i < 20; i++) {
+        auto foo = shortestPathsGenerator.getKShortest(exampleState, i);
+    }
+
     // TODO: actually write tests here
 
     EXPECT_TRUE(false);

src/utility/shortestPaths.cpp

@@ -1,7 +1,7 @@
 #include <queue>
+#include <set>
 #include "shortestPaths.h"
 #include "graph.h"
-#include "constants.h"
 
 namespace storm {
     namespace utility {
@@ -21,6 +21,8 @@ namespace storm {
 
                 // constructs the recursive shortest path representations
                 initializeShortestPaths();
+
+                candidatePaths.resize(numStates);
             }
 
             template <typename T>
@@ -53,7 +55,7 @@ namespace storm {
 
                 // set serves as priority queue with unique membership
                 // default comparison on pair actually works fine if distance is the first entry
-                std::set<std::pair<T, state_t>> dijkstraQueue;
+                std::set<std::pair<T, state_t>, std::greater_equal<std::pair<T, state_t>>> dijkstraQueue;
 
                 for (state_t initialState : model->getInitialStates()) {
                     shortestPathDistances[initialState] = zeroDistance;
@@ -124,7 +126,93 @@ namespace storm {
             }
 
             template <typename T>
-            void ShortestPathsGenerator<T>::printKShortestPath(state_t targetNode, int k, bool head) {
+            void ShortestPathsGenerator<T>::computeNextPath(state_t node, unsigned long k) {
+                assert(kShortestPaths[node].size() == k - 1); // if not, the previous SP must not exist
+
+                if (k == 2) {
+                    Path<T> shortestPathToNode = kShortestPaths[node][1 - 1]; // never forget index shift :-|
+
+                    for (state_t predecessor : graphPredecessors[node]) {
+                        // add shortest paths to predecessors plus edge to current node
+                        Path<T> pathToPredecessorPlusEdge = {
+                            boost::optional<state_t>(predecessor),
+                            1,
+                            shortestPathDistances[predecessor] * getEdgeDistance(predecessor, node)
+                        };
+                        candidatePaths[node].insert(pathToPredecessorPlusEdge);
+
+                        // ... but not the actual shortest path
+                        auto it = find(candidatePaths[node].begin(), candidatePaths[node].end(), shortestPathToNode);
+                        if (it != candidatePaths[node].end()) {
+                            candidatePaths[node].erase(it);
+                        }
+                    }
+                }
+
+                if (k > 2 || !isInitialState(node)) {
+                    // the (k-1)th shortest path (i.e., one better than the one we want to compute)
+                    Path<T> previousShortestPath = kShortestPaths[node][k - 1 - 1]; // oh god, I forgot index shift AGAIN
+
+                    // the predecessor node on that path
+                    state_t predecessor = previousShortestPath.predecessorNode.get();
+                    // the path to that predecessor was the `tailK`-shortest
+                    unsigned long tailK = previousShortestPath.predecessorK;
+
+                    // i.e. source ~~tailK-shortest path~~> predecessor --> node
+
+                    // compute one-worse-shortest path to the predecessor (if it hasn't yet been computed)
+                    if (kShortestPaths[predecessor].size() < tailK + 1) {
+                        // TODO: investigate recursion depth and possible iterative alternative
+                        computeNextPath(predecessor, tailK + 1);
+                    }
+
+                    if (kShortestPaths[predecessor].size() >= tailK + 1) {
+                        // take that path, add an edge to the current node; that's a candidate
+                        Path<T> pathToPredecessorPlusEdge = {
+                                boost::optional<state_t>(predecessor),
+                                tailK + 1,
+                                kShortestPaths[predecessor][tailK + 1 - 1].distance * getEdgeDistance(predecessor, node)
+                        };
+                        candidatePaths[node].insert(pathToPredecessorPlusEdge);
+                    }
+                    // else there was no path; TODO: does this need handling?
+                }
+
+                if (!candidatePaths[node].empty()) {
+                    Path<T> minDistanceCandidate = *(candidatePaths[node].begin());
+                    for (auto path : candidatePaths[node]) {
+                        if (path.distance > minDistanceCandidate.distance) {
+                            minDistanceCandidate = path;
+                        }
+                    }
+
+                    candidatePaths[node].erase(find(candidatePaths[node].begin(), candidatePaths[node].end(), minDistanceCandidate));
+                    kShortestPaths[node].push_back(minDistanceCandidate);
+                }
+            }
+
+            template <typename T>
+            Path<T> ShortestPathsGenerator<T>::getKShortest(state_t node, unsigned long k) {
+                unsigned long alreadyComputedK = kShortestPaths[node].size();
+
+                for (unsigned long nextK = alreadyComputedK + 1; nextK <= k; nextK++) {
+                    computeNextPath(node, nextK);
+                    if (kShortestPaths[node].size() < nextK) {
+                        break;
+                    }
+                }
+
+                if (kShortestPaths[node].size() >= k) {
+                    printKShortestPath(node, k); // DEBUG
+                } else {
+                    std::cout << "No other path exists!" << std::endl;
+                }
+
+                return kShortestPaths[node][k - 1];
+            }
+
+            template <typename T>
+            void ShortestPathsGenerator<T>::printKShortestPath(state_t targetNode, unsigned long k, bool head) {
                 // note the index shift! risk of off-by-one
                 Path<T> p = kShortestPaths[targetNode][k - 1];
 

src/utility/shortestPaths.h

@@ -5,6 +5,7 @@
 #include <boost/optional/optional.hpp>
 #include "src/models/sparse/Model.h"
 #include "src/storage/sparse/StateType.h"
+#include "constants.h"
 
 namespace storm {
     namespace utility {
@@ -26,23 +27,44 @@ namespace storm {
              * k-shortest paths list associated with t.
              *
              * Thus we can reconstruct the entire path by recursively looking
-             * up the path's tail stored as the k-th entry of the predecessor's
-             * shortest paths list.
+             * up the path's tail stored as the k-th entry [1] of the
+             * predecessor's shortest paths list.
+             *
+             * [1] oh, actually, the `k-1`th entry due to 0-based indexing!
              */
             template <typename T>
             struct Path {
                 boost::optional<state_t> predecessorNode;
-                unsigned int predecessorK;
+                unsigned long predecessorK;
                 T distance;
+
+                // FIXME: uhh.. is this okay for set? just some arbitrary order
+                bool operator<(const Path<T>& rhs) const {
+                    if (predecessorNode != rhs.predecessorNode) {
+                        return predecessorNode < rhs.predecessorNode;
+                    }
+                    return predecessorK < predecessorK;
+                }
+
+                bool operator==(const Path<T>& rhs) const {
+                    return (predecessorNode == rhs.predecessorNode) && (predecessorK == rhs.predecessorK);
+                }
             };
 
+            // -------------------------------------------------------------------------------------------------------
+
             template <typename T>
             class ShortestPathsGenerator {
             public:
                 // FIXME: this shared_ptr-passing business is probably a bad idea
                 ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model);
+
                 ~ShortestPathsGenerator();
 
+                // TODO: think about suitable output format
+                Path<T> getKShortest(state_t node, unsigned long k);
+
+
             private:
                 std::shared_ptr<storm::models::sparse::Model<T>> model;
                 storm::storage::SparseMatrix<T> transitionMatrix;
@@ -85,10 +107,32 @@ namespace storm {
                  */
                 void initializeShortestPaths();
 
+                /*!
+                 * Main step of REA algorithm. TODO: Document further.
+                 */
+                void computeNextPath(state_t node, unsigned long k);
+
                 /*!
                  * Recurses over the path and prints the nodes. Intended for debugging.
                  */
-                void printKShortestPath(state_t targetNode, int k, bool head=true);
+                void printKShortestPath(state_t targetNode, unsigned long k, bool head=true);
+
+
+                // --- tiny helper fcts ---
+                bool isInitialState(state_t node) {
+                    auto initialStates = model->getInitialStates();
+                    return find(initialStates.begin(), initialStates.end(), node) != initialStates.end();
+                }
+
+                T getEdgeDistance(state_t tailNode, state_t headNode) {
+                    for (auto const& transition : transitionMatrix.getRowGroup(tailNode)) {
+                        if (transition.getColumn() == headNode) {
+                            return transition.getValue();
+                        }
+                    }
+                    return storm::utility::zero<T>();
+                }
+                // -----------------------
             };
         }
     }