diff --git a/src/test/utility/GraphTest.cpp b/src/test/utility/GraphTest.cpp index 026f850e4..76e17299c 100644 --- a/src/test/utility/GraphTest.cpp +++ b/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> model = storm::builder::ExplicitPrismModelBuilder().translateProgram(program); - - ASSERT_TRUE(model->getType() == storm::models::ModelType::Dtmc); - - storm::storage::sparse::state_type testState = 300; - - storm::utility::shortestPaths::ShortestPathsGenerator 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); -} \ No newline at end of file diff --git a/src/utility/shortestPaths.cpp b/src/utility/shortestPaths.cpp index 8e5aef6b5..1b60717d7 100644 --- a/src/utility/shortestPaths.cpp +++ b/src/utility/shortestPaths.cpp @@ -5,12 +5,16 @@ namespace storm { namespace utility { - namespace shortestPaths { + namespace ksp { template - ShortestPathsGenerator::ShortestPathsGenerator(std::shared_ptr> model) : model(model) { + ShortestPathsGenerator::ShortestPathsGenerator(std::shared_ptr> 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(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 + ShortestPathsGenerator::ShortestPathsGenerator(std::shared_ptr> model, + state_t singleTarget) : ShortestPathsGenerator(model, std::vector{singleTarget}) {} + + template + ShortestPathsGenerator::ShortestPathsGenerator(std::shared_ptr> model, + storage::BitVector const& targetBV) : ShortestPathsGenerator(model, bitvectorToList(targetBV)) {} + template - ShortestPathsGenerator::~ShortestPathsGenerator() { + ShortestPathsGenerator::ShortestPathsGenerator(std::shared_ptr> model, + std::string const& targetLabel) : ShortestPathsGenerator(model, bitvectorToList(model->getStates(targetLabel))) {} + + template + T ShortestPathsGenerator::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(); // 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 void ShortestPathsGenerator::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 @@ -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 zeroDistance = storm::utility::one(); + T inftyDistance = utility::zero(); + T zeroDistance = utility::one(); shortestPathDistances.resize(numStates, inftyDistance); shortestPathPredecessors.resize(numStates, boost::optional()); @@ -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(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(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(currentNode); + } + // no need to enqueue meta-target } } } @@ -125,6 +179,30 @@ namespace storm { } } + template + T ShortestPathsGenerator::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(); + } else { + // edge must be "virtual edge" from target state to meta-target + assert(targetSet.count(tailNode) == 1); + return utility::one(); + } + } + + template void ShortestPathsGenerator::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 - T ShortestPathsGenerator::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(); // 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 - void ShortestPathsGenerator::printKShortestPath(state_t targetNode, unsigned long k, bool head) { + void ShortestPathsGenerator::printKShortestPath(state_t targetNode, unsigned long k, bool head) const { // note the index shift! risk of off-by-one Path 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; diff --git a/src/utility/shortestPaths.h b/src/utility/shortestPaths.h index 546afaec4..742cc2eec 100644 --- a/src/utility/shortestPaths.h +++ b/src/utility/shortestPaths.h @@ -3,6 +3,7 @@ #include #include +#include #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_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 struct Path { boost::optional 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& rhs) const { if (predecessorNode != rhs.predecessorNode) { return predecessorNode < rhs.predecessorNode; @@ -69,21 +51,36 @@ namespace storm { template class ShortestPathsGenerator { public: - // FIXME: this shared_ptr-passing business is probably a bad idea - ShortestPathsGenerator(std::shared_ptr> 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> 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> model, state_t singleTarget); + ShortestPathsGenerator(std::shared_ptr> model, storage::BitVector const& targetBV); + ShortestPathsGenerator(std::shared_ptr> 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> model; - storm::storage::SparseMatrix transitionMatrix; - state_t numStates; + std::shared_ptr> model; + storage::SparseMatrix transitionMatrix; + state_t numStates; // includes meta-target, i.e. states in model + 1 + + state_list_t targets; + std::unordered_set targetSet; + state_t metaTarget; std::vector graphPredecessors; std::vector> 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(); + return list; } // ----------------------- }; diff --git a/src/utility/shortestPaths.md b/src/utility/shortestPaths.md index 6f1042ae4..61f6fe1e3 100644 --- a/src/utility/shortestPaths.md +++ b/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. + + +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. diff --git a/test/functional/utility/KSPTest.cpp b/test/functional/utility/KSPTest.cpp new file mode 100644 index 000000000..b3118cec6 --- /dev/null +++ b/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> model = storm::builder::ExplicitPrismModelBuilder().translateProgram(program); + + storm::storage::sparse::state_type testState = 300; + storm::utility::ksp::ShortestPathsGenerator 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> model = storm::builder::ExplicitPrismModelBuilder().translateProgram(program); + + storm::storage::sparse::state_type testState = 300; + storm::utility::ksp::ShortestPathsGenerator 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> model = storm::builder::ExplicitPrismModelBuilder().translateProgram(program); + + storm::storage::sparse::state_type testState = 300; + storm::utility::ksp::ShortestPathsGenerator 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> model = storm::builder::ExplicitPrismModelBuilder().translateProgram(program); + + auto spg = storm::utility::ksp::ShortestPathsGenerator(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); +}