group targets & minimal paths

Note that this is a major, API-breaking change.

Also bunched into this commit:
 - rename namespace `shortestPaths` to `ksp`
 - omit unneeded namespace qualifiers
 - move tests from `GraphTest` to `KSPTest` and wrote more
 - path representation explanation in .md file

Former-commit-id: f395c3df40 [formerly 7fa07d1c8f]
Former-commit-id: f50dd3ce7c

src/test/utility/GraphTest.cpp

@@ -264,30 +264,3 @@ TEST(GraphTest, ExplicitProb01MinMax) {
     EXPECT_EQ(993ul, statesWithProbability01.first.getNumberOfSetBits());
     EXPECT_EQ(16ul, statesWithProbability01.second.getNumberOfSetBits());
 }
-
-
-TEST(GraphTest, kshortest) {
-    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/brp-16-2.pm");
-    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
-
-    ASSERT_TRUE(model->getType() == storm::models::ModelType::Dtmc);
-
-    storm::storage::sparse::state_type testState = 300;
-
-    storm::utility::shortestPaths::ShortestPathsGenerator<double> shortestPathsGenerator(model);
-
-    // the 1-shortest path is computed as a preprocessing step via Dijkstra;
-    // since there were some bugs here in the past, let's test it separately
-    double dijkstraSPDistance = shortestPathsGenerator.getKShortest(testState, 1);
-    std::cout << "Res: " << dijkstraSPDistance << std::endl;
-    //EXPECT_NEAR(0.0158593, dijkstraSPDistance, 0.0000001);
-    EXPECT_DOUBLE_EQ(0.015859334652581887, dijkstraSPDistance);
-
-    // main test
-    double kSPDistance1 = shortestPathsGenerator.getKShortest(testState, 100);
-    EXPECT_DOUBLE_EQ(1.5231305000339649e-06, kSPDistance1);
-
-    // let's test again to ensure re-entry is no problem
-    double kSPDistance2 = shortestPathsGenerator.getKShortest(testState, 500);
-    EXPECT_DOUBLE_EQ(3.0462610000679282e-08, kSPDistance2);
-}

src/utility/shortestPaths.cpp

@@ -5,12 +5,16 @@
 
 namespace storm {
     namespace utility {
-        namespace shortestPaths {
+        namespace ksp {
             template <typename T>
-            ShortestPathsGenerator<T>::ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model) : model(model) {
+            ShortestPathsGenerator<T>::ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model,
+                                                              state_list_t const& targets) : model(model), targets(targets) {
                 // FIXME: does this create a copy? I don't need one, so I should avoid that
                 transitionMatrix = model->getTransitionMatrix();
-                numStates = model->getNumberOfStates();
+                numStates = model->getNumberOfStates() + 1; // one more for meta-target
+
+                metaTarget = numStates - 1; // first unused state number
+                targetSet = std::unordered_set<state_t>(targets.begin(), targets.end());
 
                 computePredecessors();
 
@@ -25,22 +29,60 @@ namespace storm {
                 candidatePaths.resize(numStates);
             }
 
+            // several alternative ways to specify the targets are provided,
+            // each converts the targets and delegates to the ctor above
+            // I admit it's kind of ugly, but actually pretty convenient (and I've wanted to try C++11 delegation)
+            template <typename T>
+            ShortestPathsGenerator<T>::ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model,
+                    state_t singleTarget) : ShortestPathsGenerator<T>(model, std::vector<state_t>{singleTarget}) {}
+
+            template <typename T>
+            ShortestPathsGenerator<T>::ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model,
+                    storage::BitVector const& targetBV) : ShortestPathsGenerator<T>(model, bitvectorToList(targetBV)) {}
+
             template <typename T>
-            ShortestPathsGenerator<T>::~ShortestPathsGenerator() {
+            ShortestPathsGenerator<T>::ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model,
+                    std::string const& targetLabel) : ShortestPathsGenerator<T>(model, bitvectorToList(model->getStates(targetLabel))) {}
+
+            template <typename T>
+            T ShortestPathsGenerator<T>::computeKSP(unsigned long k) {
+                unsigned long alreadyComputedK = kShortestPaths[metaTarget].size();
+
+                for (unsigned long nextK = alreadyComputedK + 1; nextK <= k; nextK++) {
+                    computeNextPath(metaTarget, nextK);
+                    if (kShortestPaths[metaTarget].size() < nextK) {
+                        std::cout << std::endl << "--> DEBUG: Last path: k=" << (nextK - 1) << ":" << std::endl;
+                        printKShortestPath(metaTarget, nextK - 1);
+                        std::cout << "---------" << "No other path exists!" << std::endl;
+                        return utility::zero<T>(); // TODO: throw exception or something
+                    }
+                }
+
+                //std::cout << std::endl << "--> DEBUG: Finished. " << k << "-shortest path:" << std::endl;
+                //printKShortestPath(metaTarget, k);
+                //std::cout << "---------" << std::endl;
+                return kShortestPaths[metaTarget][k - 1].distance;
             }
 
             template <typename T>
             void ShortestPathsGenerator<T>::computePredecessors() {
-                auto rowCount = transitionMatrix.getRowCount();
-                assert(numStates == rowCount);
+                assert(numStates - 1 == transitionMatrix.getRowCount());
 
-                graphPredecessors.resize(rowCount);
+                // one more for meta-target
+                graphPredecessors.resize(numStates);
 
-                for (state_t i = 0; i < rowCount; i++) {
+                for (state_t i = 0; i < numStates - 1; i++) {
                     for (auto const& transition : transitionMatrix.getRowGroup(i)) {
-                        graphPredecessors[transition.getColumn()].push_back(i);
+                        // to avoid non-minimal paths, the target states are
+                        // *not* predecessors of any state but the meta-target
+                        if (std::find(targets.begin(), targets.end(), i) == targets.end()) {
+                            graphPredecessors[transition.getColumn()].push_back(i);
+                        }
                     }
                 }
+
+                // meta-target has exactly the target states as predecessors
+                graphPredecessors[metaTarget] = targets;
             }
 
             template <typename T>
@@ -48,8 +90,8 @@ namespace storm {
                 // the existing Dijkstra isn't working anyway AND
                 // doesn't fully meet our requirements, so let's roll our own
 
-                T inftyDistance = storm::utility::zero<T>();
-                T zeroDistance = storm::utility::one<T>();
+                T inftyDistance = utility::zero<T>();
+                T zeroDistance = utility::one<T>();
                 shortestPathDistances.resize(numStates, inftyDistance);
                 shortestPathPredecessors.resize(numStates, boost::optional<state_t>());
 
@@ -66,16 +108,28 @@ namespace storm {
                     state_t currentNode = (*dijkstraQueue.begin()).second;
                     dijkstraQueue.erase(dijkstraQueue.begin());
 
-                    for (auto const& transition : transitionMatrix.getRowGroup(currentNode)) {
-                        state_t otherNode = transition.getColumn();
-
-                        // note that distances are probabilities, thus they are multiplied and larger is better
-                        T alternateDistance = shortestPathDistances[currentNode] * transition.getValue();
-                        if (alternateDistance > shortestPathDistances[otherNode]) {
-                            shortestPathDistances[otherNode] = alternateDistance;
-                            shortestPathPredecessors[otherNode] = boost::optional<state_t>(currentNode);
-                            dijkstraQueue.emplace(alternateDistance, otherNode);
+                    if (targetSet.count(currentNode) == 0) {
+                        // non-target node, treated normally
+                        for (auto const& transition : transitionMatrix.getRowGroup(currentNode)) {
+                            state_t otherNode = transition.getColumn();
+
+                            // note that distances are probabilities, thus they are multiplied and larger is better
+                            T alternateDistance = shortestPathDistances[currentNode] * transition.getValue();
+                            if (alternateDistance > shortestPathDistances[otherNode]) {
+                                shortestPathDistances[otherNode] = alternateDistance;
+                                shortestPathPredecessors[otherNode] = boost::optional<state_t>(currentNode);
+                                dijkstraQueue.emplace(alternateDistance, otherNode);
+                            }
                         }
+                    } else {
+                        // target node has only "virtual edge" (with prob 1) to meta-target
+                        // no multiplication necessary
+                        T alternateDistance = shortestPathDistances[currentNode];
+                        if (alternateDistance > shortestPathDistances[metaTarget]) {
+                            shortestPathDistances[metaTarget] = alternateDistance;
+                            shortestPathPredecessors[metaTarget] = boost::optional<state_t>(currentNode);
+                        }
+                        // no need to enqueue meta-target
                     }
                 }
             }
@@ -125,6 +179,30 @@ namespace storm {
                 }
             }
 
+            template <typename T>
+            T ShortestPathsGenerator<T>::getEdgeDistance(state_t tailNode, state_t headNode) const {
+                // just to be clear, head is where the arrow points (obviously)
+                if (headNode != metaTarget) {
+                    // edge is "normal", not to meta-target
+
+                    for (auto const& transition : transitionMatrix.getRowGroup(tailNode)) {
+                        if (transition.getColumn() == headNode) {
+                            return transition.getValue();
+                        }
+                    }
+
+                    // there is no such edge
+                    // let's disallow that for now, because I'm not expecting it to happen
+                    assert(false);
+                    return utility::zero<T>();
+                } else {
+                    // edge must be "virtual edge" from target state to meta-target
+                    assert(targetSet.count(tailNode) == 1);
+                    return utility::one<T>();
+                }
+            }
+
+
             template <typename T>
             void ShortestPathsGenerator<T>::computeNextPath(state_t node, unsigned long k) {
                 assert(kShortestPaths[node].size() == k - 1); // if not, the previous SP must not exist
@@ -142,7 +220,7 @@ namespace storm {
                         candidatePaths[node].insert(pathToPredecessorPlusEdge);
 
                         // ... but not the actual shortest path
-                        auto it = find(candidatePaths[node].begin(), candidatePaths[node].end(), shortestPathToNode);
+                        auto it = std::find(candidatePaths[node].begin(), candidatePaths[node].end(), shortestPathToNode);
                         if (it != candidatePaths[node].end()) {
                             candidatePaths[node].erase(it);
                         }
@@ -186,38 +264,22 @@ namespace storm {
                         }
                     }
 
-                    candidatePaths[node].erase(find(candidatePaths[node].begin(), candidatePaths[node].end(), minDistanceCandidate));
+                    candidatePaths[node].erase(std::find(candidatePaths[node].begin(), candidatePaths[node].end(), minDistanceCandidate));
                     kShortestPaths[node].push_back(minDistanceCandidate);
                 }
             }
 
             template <typename T>
-            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) {
-                        std::cout << std::endl << "--> DEBUG: Last path: k=" << (nextK - 1) << ":" << std::endl;
-                        printKShortestPath(node, nextK - 1);
-                        std::cout << "---------" << "No other path exists!" << std::endl;
-                        return storm::utility::zero<T>(); // TODO: throw exception or something
-                    }
-                }
-
-                std::cout << std::endl << "--> DEBUG: Finished. " << k << "-shortest path:" << std::endl;
-                printKShortestPath(node, k);
-                std::cout << "---------" << std::endl;
-                return kShortestPaths[node][k - 1].distance;
-            }
-
-            template <typename T>
-            void ShortestPathsGenerator<T>::printKShortestPath(state_t targetNode, unsigned long k, bool head) {
+            void ShortestPathsGenerator<T>::printKShortestPath(state_t targetNode, unsigned long k, bool head) const {
                 // note the index shift! risk of off-by-one
                 Path<T> p = kShortestPaths[targetNode][k - 1];
 
                 if (head) {
-                    std::cout << "Path (reversed), dist (prob)=" << p.distance << ": [";
+                    std::cout << "Path (reversed";
+                    if (targetNode == metaTarget) {
+                        std::cout << ", w/ meta-target";
+                    }
+                    std::cout <<"), dist (prob)=" << p.distance << ": [";
                 }
 
                 std::cout << " " << targetNode;

src/utility/shortestPaths.h

@@ -3,6 +3,7 @@
 
 #include <vector>
 #include <boost/optional/optional.hpp>
+#include <unordered_set>
 #include "src/models/sparse/Model.h"
 #include "src/storage/sparse/StateType.h"
 #include "constants.h"
@@ -22,36 +23,17 @@
 
 namespace storm {
     namespace utility {
-        namespace shortestPaths {
-            typedef storm::storage::sparse::state_type state_t;
+        namespace ksp {
+            typedef storage::sparse::state_type state_t;
             typedef std::vector<state_t> state_list_t;
 
-            /*
-             * Implicit shortest path representation
-             *
-             * All shortest paths (from s to t) can be described as some
-             * k-shortest path to some node u plus the edge to t:
-             *
-             *     s ~~k-shortest path~~> u --> t
-             *
-             * This struct stores u (`predecessorNode`) and k (`predecessorK`).
-             *
-             * t is implied by this struct's location: It is stored in the
-             * 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 [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 long predecessorK;
                 T distance;
 
-                // FIXME: uhh.. is this okay for set? just some arbitrary order
+                // arbitrary order for std::set
                 bool operator<(const Path<T>& rhs) const {
                     if (predecessorNode != rhs.predecessorNode) {
                         return predecessorNode < rhs.predecessorNode;
@@ -69,21 +51,36 @@ namespace storm {
             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);
+                /*!
+                 * Performs precomputations (including meta-target insertion and Dijkstra).
+                 * Modifications are done locally, `model` remains unchanged.
+                 * Target (group) cannot be changed.
+                 */
+                // FIXME: this shared_ptr-passing business might be a bad idea
+                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, state_list_t const& targets);
 
-                ~ShortestPathsGenerator();
+                // allow alternative ways of specifying the target,
+                // all of which will be converted to list and delegated to constructor above
+                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, state_t singleTarget);
+                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, storage::BitVector const& targetBV);
+                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, std::string const& targetLabel);
 
-                // TODO: think about suitable output format
-                T getKShortest(state_t node, unsigned long k);
+                ~ShortestPathsGenerator(){}
 
-                // TODO: allow three kinds of target arguments: single state, BitVector, Label
+                /*!
+                 * Computes k-SP to target and returns distance (probability).
+                 */
+                T computeKSP(unsigned long k);
 
 
             private:
-                std::shared_ptr<storm::models::sparse::Model<T>> model;
-                storm::storage::SparseMatrix<T> transitionMatrix;
-                state_t numStates;
+                std::shared_ptr<models::sparse::Model<T>> model;
+                storage::SparseMatrix<T> transitionMatrix;
+                state_t numStates; // includes meta-target, i.e. states in model + 1
+
+                state_list_t targets;
+                std::unordered_set<state_t> targetSet;
+                state_t metaTarget;
 
                 std::vector<state_list_t>             graphPredecessors;
                 std::vector<boost::optional<state_t>> shortestPathPredecessors;
@@ -130,22 +127,26 @@ namespace storm {
                 /*!
                  * Recurses over the path and prints the nodes. Intended for debugging.
                  */
-                void printKShortestPath(state_t targetNode, unsigned long k, bool head=true);
+                void printKShortestPath(state_t targetNode, unsigned long k, bool head=true) const;
+
+                /*!
+                 * Returns actual distance for real edges, 1 for edges to meta-target.
+                 */
+                T getEdgeDistance(state_t tailNode, state_t headNode) const;
 
 
                 // --- tiny helper fcts ---
-                bool isInitialState(state_t node) {
+                bool isInitialState(state_t node) const {
                     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();
-                        }
+                state_list_t bitvectorToList(storage::BitVector const& bv) const {
+                    state_list_t list;
+                    for (state_t state : bv) {
+                        list.push_back(state);
                     }
-                    return storm::utility::zero<T>();
+                    return list;
                 }
                 // -----------------------
             };

src/utility/shortestPaths.md

@@ -4,7 +4,8 @@ the design decisions and the rationale behind them.]
 
 ## Differences from REA algorithm
 This class closely follows the Jimenez-Marzal REA algorithm.
-However, there are some notable deviations:
+However, there are some notable deviations in the way targets and shortest
+paths are defined.
 
 ### Target groups
 Firstly, instead of a single target state, a group of target states is
@@ -13,6 +14,7 @@ transitions (other than self-loops, but this is moot due to not allowing
 non-minimal paths -- see below) and instead introducing edges to a new sink
 state that serves as a meta-target.
 
+<!--
 In terms of implementation, there are two possible routes (that I can think
 of):
 
@@ -29,11 +31,18 @@ of):
 It is not clear if there will ever be a need for changing targets. While
 the overlay option is alluring, in the spirit of YAGNI, I choose the
 destructive, simple route.
+-->
+
+I chose to implement this by modifying the precomputation results, meaning
+that they are only valid for a fixed group of target states. Thus, the
+target group is required in the constructor. (It would have been possible to
+allow for interchangeable target groups, but I don't anticipate that use
+case.)
 
 ### Minimality of paths
 Secondly, we define shortest paths as "minimal" shortest paths in the
-following sense: The path may only visit the target state once (at the
-end). In other words, no KSP may be a prefix of another.
+following sense: The path may not visit any target state except at the
+end. As a corollary, no KSP (to a target node) may be a prefix of another.
 This in particular forbids shortest path progressions such as these:
 
     1-SP: 1 2 3
@@ -48,3 +57,30 @@ interested in such paths.
 shorter paths is non-empty (which is an even stronger statement than
 loop-free-ness of paths), because we want to take a union of those node
 sets. But that's a different matter.)
+
+
+## Data structures, in particular: Path representation
+
+The implicit shortest path representation that J&M describe in the paper
+is used, except that indices rather than back-pointers are stored to
+refer to the tail of the path.
+[Maybe pointers would be faster? STL vector access via index should be
+pretty fast too, though, and less error-prone.]
+
+A bit more detail (recap of the paper):
+All shortest paths (from `s` to `t`) can be described as some k-shortest
+path to some node `u` plus an edge to `t`:
+
+    s ~~k-shortest path~~> u --> t
+
+Further, the shortest paths to some node are always computed in order and
+without gaps, e.g., the 1, 2, 3-shortest paths to `t` will be computed
+before the 4-SP. Thus, we store the SPs in a linked list for each node,
+with the k-th entry[^1] being the k-th SP to that node.
+
+Thus for an SP as shown above we simply store the predecessor node (`u`)
+and the `k`, which allows us to look up the tail of the SP.
+By recursively looking up the tail (until it's empty), we reconstruct
+the entire path back-to-front.
+
+[^1]: Which due to 0-based indexing has index `k-1`, of course! Damn it.

test/functional/utility/KSPTest.cpp

@@ -0,0 +1,69 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/parser/PrismParser.h"
+#include "src/models/sparse/Dtmc.h"
+#include "src/builder/ExplicitPrismModelBuilder.h"
+#include "src/utility/graph.h"
+#include "src/utility/shortestPaths.cpp"
+
+// NOTE: These result values of these tests were generated by the
+//       KSP-Generator itself and checked for gross implausibility, but no
+//       more than that.
+//       An independent verification of the values would be really nice ...
+
+TEST(KSPTest, dijkstra) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/brp-16-2.pm");
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    storm::storage::sparse::state_type testState = 300;
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+
+    double dist = spg.computeKSP(1);
+    EXPECT_DOUBLE_EQ(0.015859334652581887, dist);
+}
+
+TEST(KSPTest, singleTarget) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/brp-16-2.pm");
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    storm::storage::sparse::state_type testState = 300;
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+
+    double dist = spg.computeKSP(100);
+    EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist);
+}
+
+TEST(KSPTest, reentry) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/brp-16-2.pm");
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    storm::storage::sparse::state_type testState = 300;
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+
+    double dist = spg.computeKSP(100);
+    EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist);
+
+    // get another distance to ensure re-entry is no problem
+    double dist2 = spg.computeKSP(500);
+    EXPECT_DOUBLE_EQ(3.0462610000679282e-08, dist2);
+}
+
+TEST(KSPTest, groupTarget) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/brp-16-2.pm");
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>().translateProgram(program);
+
+    auto spg = storm::utility::ksp::ShortestPathsGenerator<double>(model, storm::utility::ksp::state_list_t{50, 90});
+
+    // this path should lead to 90
+    double dist1 = spg.computeKSP(8);
+    EXPECT_DOUBLE_EQ(7.3796982335999988e-06, dist1);
+
+    // this one to 50
+    double dist2 = spg.computeKSP(9);
+    EXPECT_DOUBLE_EQ(3.92e-06, dist2);
+
+    // this one to 90 again
+    double dist3 = spg.computeKSP(12);
+    EXPECT_DOUBLE_EQ(3.6160521344640002e-06, dist3);
+}