dehnert
8 years ago
9 changed files with 741 additions and 92 deletions
-
4src/storm/utility/eigen.h
-
1src/storm/utility/gmm.h
-
69src/storm/utility/graph.cpp
-
29src/storm/utility/graph.h
-
376src/storm/utility/shortestPaths.cpp
-
243src/storm/utility/shortestPaths.h
-
1src/test/utility/GraphTest.cpp
-
108src/test/utility/KSPTest.cpp
-
2storm-config.h.in
@ -0,0 +1,376 @@ |
|||
#include <queue>
|
|||
#include <set>
|
|||
#include <string>
|
|||
|
|||
#include "macros.h"
|
|||
#include "shortestPaths.h"
|
|||
#include "graph.h"
|
|||
|
|||
// FIXME: I've accidentally used k=0 *twice* now without realizing that k>=1 is required!
|
|||
// Accessing zero should trigger a warning!
|
|||
// (Also, did I document this? I think so, somewhere. I went with k>=1 because
|
|||
// that's what the KSP paper used, but in retrospect k>=0 seems more intuitive ...)
|
|||
|
|||
namespace storm { |
|||
namespace utility { |
|||
namespace ksp { |
|||
template <typename T> |
|||
ShortestPathsGenerator<T>::ShortestPathsGenerator(storage::SparseMatrix<T> const& transitionMatrix, |
|||
std::unordered_map<state_t, T> const& targetProbMap, |
|||
BitVector const& initialStates, |
|||
MatrixFormat matrixFormat) : |
|||
transitionMatrix(transitionMatrix), |
|||
numStates(transitionMatrix.getColumnCount() + 1), // one more for meta-target
|
|||
metaTarget(transitionMatrix.getColumnCount()), // first unused state index
|
|||
initialStates(initialStates), |
|||
targetProbMap(targetProbMap), |
|||
matrixFormat(matrixFormat) { |
|||
|
|||
computePredecessors(); |
|||
|
|||
// gives us SP-predecessors, SP-distances
|
|||
performDijkstra(); |
|||
|
|||
computeSPSuccessors(); |
|||
|
|||
// constructs the recursive shortest path representations
|
|||
initializeShortestPaths(); |
|||
|
|||
candidatePaths.resize(numStates); |
|||
} |
|||
|
|||
template <typename T> |
|||
ShortestPathsGenerator<T>::ShortestPathsGenerator(storage::SparseMatrix<T> const& transitionMatrix, |
|||
std::vector<T> const& targetProbVector, BitVector const& initialStates, MatrixFormat matrixFormat) |
|||
: ShortestPathsGenerator<T>(transitionMatrix, vectorToMap(targetProbVector), initialStates, matrixFormat) {} |
|||
|
|||
// extracts the relevant info from the model and delegates to ctor above
|
|||
template <typename T> |
|||
ShortestPathsGenerator<T>::ShortestPathsGenerator(Model const& model, BitVector const& targetBV) |
|||
: ShortestPathsGenerator<T>(model.getTransitionMatrix(), allProbOneMap(targetBV), model.getInitialStates(), MatrixFormat::straight) {} |
|||
|
|||
// 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(Model const& model, state_t singleTarget) |
|||
: ShortestPathsGenerator<T>(model, std::vector<state_t>{singleTarget}) {} |
|||
|
|||
template <typename T> |
|||
ShortestPathsGenerator<T>::ShortestPathsGenerator(Model const& model, std::vector<state_t> const& targetList) |
|||
: ShortestPathsGenerator<T>(model, BitVector(model.getNumberOfStates(), targetList)) {} |
|||
|
|||
template <typename T> |
|||
ShortestPathsGenerator<T>::ShortestPathsGenerator(Model const& model, std::string const& targetLabel) |
|||
: ShortestPathsGenerator<T>(model, model.getStates(targetLabel)) {} |
|||
|
|||
template <typename T> |
|||
T ShortestPathsGenerator<T>::getDistance(unsigned long k) { |
|||
computeKSP(k); |
|||
return kShortestPaths[metaTarget][k - 1].distance; |
|||
} |
|||
|
|||
template <typename T> |
|||
BitVector ShortestPathsGenerator<T>::getStates(unsigned long k) { |
|||
computeKSP(k); |
|||
BitVector stateSet(numStates - 1, false); // no meta-target
|
|||
|
|||
Path<T> currentPath = kShortestPaths[metaTarget][k - 1]; |
|||
boost::optional<state_t> maybePredecessor = currentPath.predecessorNode; |
|||
// this omits the first node, which is actually convenient since that's the meta-target
|
|||
|
|||
while (maybePredecessor) { |
|||
state_t predecessor = maybePredecessor.get(); |
|||
stateSet.set(predecessor, true); |
|||
|
|||
currentPath = kShortestPaths[predecessor][currentPath.predecessorK - 1]; // god damn you, index
|
|||
maybePredecessor = currentPath.predecessorNode; |
|||
} |
|||
|
|||
return stateSet; |
|||
} |
|||
|
|||
template <typename T> |
|||
std::vector<state_t> ShortestPathsGenerator<T>::getPathAsList(unsigned long k) { |
|||
computeKSP(k); |
|||
|
|||
std::vector<state_t> backToFrontList; |
|||
|
|||
Path<T> currentPath = kShortestPaths[metaTarget][k - 1]; |
|||
boost::optional<state_t> maybePredecessor = currentPath.predecessorNode; |
|||
// this omits the first node, which is actually convenient since that's the meta-target
|
|||
|
|||
while (maybePredecessor) { |
|||
state_t predecessor = maybePredecessor.get(); |
|||
backToFrontList.push_back(predecessor); |
|||
|
|||
currentPath = kShortestPaths[predecessor][currentPath.predecessorK - 1]; |
|||
maybePredecessor = currentPath.predecessorNode; |
|||
} |
|||
|
|||
return backToFrontList; |
|||
} |
|||
|
|||
template <typename T> |
|||
void ShortestPathsGenerator<T>::computePredecessors() { |
|||
assert(transitionMatrix.hasTrivialRowGrouping()); |
|||
|
|||
// one more for meta-target
|
|||
graphPredecessors.resize(numStates); |
|||
|
|||
for (state_t i = 0; i < numStates - 1; i++) { |
|||
for (auto const& transition : transitionMatrix.getRowGroup(i)) { |
|||
// to avoid non-minimal paths, the meta-target-predecessors are
|
|||
// *not* predecessors of any state but the meta-target
|
|||
if (!isMetaTargetPredecessor(i)) { |
|||
graphPredecessors[transition.getColumn()].push_back(i); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// meta-target has exactly the meta-target-predecessors as predecessors
|
|||
// (duh. note that the meta-target-predecessors used to be called target,
|
|||
// but that's not necessarily true in the matrix/value invocation case)
|
|||
for (auto targetProbPair : targetProbMap) { |
|||
graphPredecessors[metaTarget].push_back(targetProbPair.first); |
|||
} |
|||
} |
|||
|
|||
template <typename T> |
|||
void ShortestPathsGenerator<T>::performDijkstra() { |
|||
// the existing Dijkstra isn't working anyway AND
|
|||
// doesn't fully meet our requirements, so let's roll our own
|
|||
|
|||
T inftyDistance = zero<T>(); |
|||
T zeroDistance = one<T>(); |
|||
shortestPathDistances.resize(numStates, inftyDistance); |
|||
shortestPathPredecessors.resize(numStates, boost::optional<state_t>()); |
|||
|
|||
// 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>, std::greater<std::pair<T, state_t>>> dijkstraQueue; |
|||
|
|||
for (state_t initialState : initialStates) { |
|||
shortestPathDistances[initialState] = zeroDistance; |
|||
dijkstraQueue.emplace(zeroDistance, initialState); |
|||
} |
|||
|
|||
while (!dijkstraQueue.empty()) { |
|||
state_t currentNode = (*dijkstraQueue.begin()).second; |
|||
dijkstraQueue.erase(dijkstraQueue.begin()); |
|||
|
|||
if (!isMetaTargetPredecessor(currentNode)) { |
|||
// 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] * convertDistance(currentNode, otherNode, transition.getValue()); |
|||
assert(zero<T>() <= alternateDistance <= one<T>()); // FIXME: there is a negative transition! SM gives us a placeholder!
|
|||
if (alternateDistance > shortestPathDistances[otherNode]) { |
|||
shortestPathDistances[otherNode] = alternateDistance; |
|||
shortestPathPredecessors[otherNode] = boost::optional<state_t>(currentNode); |
|||
dijkstraQueue.emplace(alternateDistance, otherNode); |
|||
} |
|||
} |
|||
} else { |
|||
// node only has one "virtual edge" (with prob as per targetProbMap) to meta-target
|
|||
// FIXME: double check
|
|||
T alternateDistance = shortestPathDistances[currentNode] * targetProbMap[currentNode]; |
|||
if (alternateDistance > shortestPathDistances[metaTarget]) { |
|||
shortestPathDistances[metaTarget] = alternateDistance; |
|||
shortestPathPredecessors[metaTarget] = boost::optional<state_t>(currentNode); |
|||
} |
|||
// no need to enqueue meta-target
|
|||
} |
|||
} |
|||
} |
|||
|
|||
template <typename T> |
|||
void ShortestPathsGenerator<T>::computeSPSuccessors() { |
|||
shortestPathSuccessors.resize(numStates); |
|||
|
|||
for (state_t i = 0; i < numStates; i++) { |
|||
if (shortestPathPredecessors[i]) { |
|||
state_t predecessor = shortestPathPredecessors[i].get(); |
|||
shortestPathSuccessors[predecessor].push_back(i); |
|||
} |
|||
} |
|||
} |
|||
|
|||
template <typename T> |
|||
void ShortestPathsGenerator<T>::initializeShortestPaths() { |
|||
kShortestPaths.resize(numStates); |
|||
|
|||
// BFS in Dijkstra-SP order
|
|||
std::queue<state_t> bfsQueue; |
|||
for (state_t initialState : initialStates) { |
|||
bfsQueue.push(initialState); |
|||
} |
|||
|
|||
while (!bfsQueue.empty()) { |
|||
state_t currentNode = bfsQueue.front(); |
|||
bfsQueue.pop(); |
|||
|
|||
if (!kShortestPaths[currentNode].empty()) { |
|||
continue; // already visited
|
|||
} |
|||
|
|||
for (state_t successorNode : shortestPathSuccessors[currentNode]) { |
|||
bfsQueue.push(successorNode); |
|||
} |
|||
|
|||
// note that `shortestPathPredecessor` may not be present
|
|||
// if current node is an initial state
|
|||
// in this case, the boost::optional copy of an uninitialized optional is hopefully also uninitialized
|
|||
kShortestPaths[currentNode].push_back(Path<T> { |
|||
shortestPathPredecessors[currentNode], |
|||
1, |
|||
shortestPathDistances[currentNode] |
|||
}); |
|||
} |
|||
} |
|||
|
|||
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 convertDistance(tailNode, headNode, transition.getValue()); |
|||
} |
|||
} |
|||
|
|||
// there is no such edge
|
|||
// let's disallow that for now, because I'm not expecting it to happen
|
|||
assert(false); |
|||
return zero<T>(); |
|||
} else { |
|||
// edge must be "virtual edge" to meta-target
|
|||
assert(isMetaTargetPredecessor(tailNode)); |
|||
return targetProbMap.at(tailNode); // FIXME double check
|
|||
} |
|||
} |
|||
|
|||
|
|||
template <typename T> |
|||
void ShortestPathsGenerator<T>::computeNextPath(state_t node, unsigned long k) { |
|||
assert(k >= 2); // Dijkstra is used for k=1
|
|||
assert(kShortestPaths[node].size() == k - 1); // if not, the previous SP must not exist
|
|||
|
|||
// TODO: I could extract the candidate generation to make this function more succinct
|
|||
if (k == 2) { |
|||
// Step B.1 in J&M paper
|
|||
|
|||
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 = std::find(candidatePaths[node].begin(), candidatePaths[node].end(), shortestPathToNode); |
|||
if (it != candidatePaths[node].end()) { |
|||
candidatePaths[node].erase(it); |
|||
} |
|||
} |
|||
} |
|||
|
|||
if (not (k == 2 && isInitialState(node))) { |
|||
// Steps B.2-5 in J&M paper
|
|||
|
|||
// 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? -- yes, but not here (because the step B.1 may have added candidates)
|
|||
} |
|||
|
|||
// Step B.6 in J&M paper
|
|||
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(std::find(candidatePaths[node].begin(), candidatePaths[node].end(), minDistanceCandidate)); |
|||
kShortestPaths[node].push_back(minDistanceCandidate); |
|||
} else { |
|||
// TODO: kSP does not exist. this is handled later, but it would be nice to catch it as early as possble, wouldn't it?
|
|||
STORM_LOG_TRACE("KSP: no candidates, this will trigger nonexisting ksp after exiting these recursions. TODO: handle here"); |
|||
} |
|||
} |
|||
|
|||
template <typename T> |
|||
void 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) { |
|||
unsigned long lastExistingK = nextK - 1; |
|||
STORM_LOG_DEBUG("KSP throws (as expected) due to nonexistence -- maybe this is unhandled and causes the Python interface to segfault?"); |
|||
STORM_LOG_DEBUG("last existing k-SP has k=" + std::to_string(lastExistingK)); |
|||
STORM_LOG_DEBUG("maybe this is unhandled and causes the Python interface to segfault?"); |
|||
throw std::invalid_argument("k-SP does not exist for k=" + std::to_string(k)); |
|||
} |
|||
} |
|||
} |
|||
|
|||
template <typename T> |
|||
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"; |
|||
if (targetNode == metaTarget) { |
|||
std::cout << ", w/ meta-target"; |
|||
} |
|||
std::cout <<"), dist (prob)=" << p.distance << ": ["; |
|||
} |
|||
|
|||
std::cout << " " << targetNode; |
|||
|
|||
if (p.predecessorNode) { |
|||
printKShortestPath(p.predecessorNode.get(), p.predecessorK, false); |
|||
} else { |
|||
std::cout << " ]" << std::endl; |
|||
} |
|||
} |
|||
|
|||
|
|||
template class ShortestPathsGenerator<double>; |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,243 @@ |
|||
#ifndef STORM_UTIL_SHORTESTPATHS_H_ |
|||
#define STORM_UTIL_SHORTESTPATHS_H_ |
|||
|
|||
#include <vector> |
|||
#include <boost/optional/optional.hpp> |
|||
#include <unordered_set> |
|||
#include "storm/models/sparse/Model.h" |
|||
#include "storm/storage/sparse/StateType.h" |
|||
#include "constants.h" |
|||
|
|||
// NOTE: You'll (eventually) find the usual API documentation below; |
|||
// for more information about the purpose, design decisions, |
|||
// etc., you may consult `shortestPath.md`. - Tom |
|||
|
|||
// TODO: test whether using BitVector instead of vector<state_t> is |
|||
// faster for storing predecessors etc. |
|||
|
|||
// Now using BitVectors instead of vector<state_t> in the API because |
|||
// BitVectors are used throughout Storm to represent a (unordered) list |
|||
// of states. |
|||
// (Even though both initialStates and targets are probably very sparse.) |
|||
|
|||
namespace storm { |
|||
namespace utility { |
|||
namespace ksp { |
|||
using state_t = storage::sparse::state_type; |
|||
using OrderedStateList = std::vector<state_t>; |
|||
using BitVector = storage::BitVector; |
|||
|
|||
// -- helper structs/classes ----------------------------------------------------------------------------- |
|||
|
|||
template <typename T> |
|||
struct Path { |
|||
boost::optional<state_t> predecessorNode; |
|||
unsigned long predecessorK; |
|||
T distance; |
|||
|
|||
// arbitrary order for std::set |
|||
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); |
|||
} |
|||
}; |
|||
|
|||
// when using the raw matrix/vector invocation, this enum parameter |
|||
// forces the caller to declare whether the matrix has the evil I-P |
|||
// format, which requires back-conversion of the entries |
|||
enum class MatrixFormat { straight, iMinusP }; |
|||
|
|||
// ------------------------------------------------------------------------------------------------------- |
|||
|
|||
template <typename T> |
|||
class ShortestPathsGenerator { |
|||
public: |
|||
using Matrix = storage::SparseMatrix<T>; |
|||
using StateProbMap = std::unordered_map<state_t, T>; |
|||
using Model = models::sparse::Model<T>; |
|||
|
|||
/*! |
|||
* Performs precomputations (including meta-target insertion and Dijkstra). |
|||
* Modifications are done locally, `model` remains unchanged. |
|||
* Target (group) cannot be changed. |
|||
*/ |
|||
ShortestPathsGenerator(Model const& model, BitVector const& targetBV); |
|||
|
|||
// allow alternative ways of specifying the target, |
|||
// all of which will be converted to BitVector and delegated to constructor above |
|||
ShortestPathsGenerator(Model const& model, state_t singleTarget); |
|||
ShortestPathsGenerator(Model const& model, std::vector<state_t> const& targetList); |
|||
ShortestPathsGenerator(Model const& model, std::string const& targetLabel = "target"); |
|||
|
|||
// a further alternative: use transition matrix of maybe-states |
|||
// combined with target vector (e.g., the instantiated matrix/vector from SamplingModel); |
|||
// in this case separately specifying a target makes no sense |
|||
ShortestPathsGenerator(Matrix const& transitionMatrix, std::vector<T> const& targetProbVector, BitVector const& initialStates, MatrixFormat matrixFormat); |
|||
ShortestPathsGenerator(Matrix const& maybeTransitionMatrix, StateProbMap const& targetProbMap, BitVector const& initialStates, MatrixFormat matrixFormat); |
|||
|
|||
|
|||
inline ~ShortestPathsGenerator(){} |
|||
|
|||
/*! |
|||
* Returns distance (i.e., probability) of the KSP. |
|||
* Computes KSP if not yet computed. |
|||
* @throws std::invalid_argument if no such k-shortest path exists |
|||
*/ |
|||
T getDistance(unsigned long k); |
|||
|
|||
/*! |
|||
* Returns the states that occur in the KSP. |
|||
* For a path-traversal (in order and with duplicates), see `getPathAsList`. |
|||
* Computes KSP if not yet computed. |
|||
* @throws std::invalid_argument if no such k-shortest path exists |
|||
*/ |
|||
storage::BitVector getStates(unsigned long k); |
|||
|
|||
/*! |
|||
* Returns the states of the KSP as back-to-front traversal. |
|||
* Computes KSP if not yet computed. |
|||
* @throws std::invalid_argument if no such k-shortest path exists |
|||
*/ |
|||
OrderedStateList getPathAsList(unsigned long k); |
|||
|
|||
|
|||
private: |
|||
Matrix const& transitionMatrix; |
|||
state_t numStates; // includes meta-target, i.e. states in model + 1 |
|||
state_t metaTarget; |
|||
BitVector initialStates; |
|||
StateProbMap targetProbMap; |
|||
|
|||
MatrixFormat matrixFormat; |
|||
|
|||
std::vector<OrderedStateList> graphPredecessors; |
|||
std::vector<boost::optional<state_t>> shortestPathPredecessors; |
|||
std::vector<OrderedStateList> shortestPathSuccessors; |
|||
std::vector<T> shortestPathDistances; |
|||
|
|||
std::vector<std::vector<Path<T>>> kShortestPaths; |
|||
std::vector<std::set<Path<T>>> candidatePaths; |
|||
|
|||
/*! |
|||
* Computes list of predecessors for all nodes. |
|||
* Reachability is not considered; a predecessor is simply any node that has an edge leading to the node in question. |
|||
* Requires `transitionMatrix`. |
|||
* Modifies `graphPredecessors`. |
|||
*/ |
|||
void computePredecessors(); |
|||
|
|||
/*! |
|||
* Computes shortest path distances and predecessors. |
|||
* Requires `model`, `numStates`, `transitionMatrix`. |
|||
* Modifies `shortestPathPredecessors` and `shortestPathDistances`. |
|||
*/ |
|||
void performDijkstra(); |
|||
|
|||
/*! |
|||
* Computes list of shortest path successor nodes from predecessor list. |
|||
* Requires `shortestPathPredecessors`, `numStates`. |
|||
* Modifies `shortestPathSuccessors`. |
|||
*/ |
|||
void computeSPSuccessors(); |
|||
|
|||
/*! |
|||
* Constructs and stores the implicit shortest path representations (see `Path`) for the (1-)shortest paths. |
|||
* Requires `shortestPathPredecessors`, `shortestPathDistances`, `model`, `numStates`. |
|||
* Modifies `kShortestPaths`. |
|||
*/ |
|||
void initializeShortestPaths(); |
|||
|
|||
/*! |
|||
* Main step of REA algorithm. TODO: Document further. |
|||
*/ |
|||
void computeNextPath(state_t node, unsigned long k); |
|||
|
|||
/*! |
|||
* Computes k-shortest path if not yet computed. |
|||
* @throws std::invalid_argument if no such k-shortest path exists |
|||
*/ |
|||
void computeKSP(unsigned long k); |
|||
|
|||
/*! |
|||
* Recurses over the path and prints the nodes. Intended for debugging. |
|||
*/ |
|||
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 --- |
|||
|
|||
inline bool isInitialState(state_t node) const { |
|||
return find(initialStates.begin(), initialStates.end(), node) != initialStates.end(); |
|||
} |
|||
|
|||
inline bool isMetaTargetPredecessor(state_t node) const { |
|||
return targetProbMap.count(node) == 1; |
|||
} |
|||
|
|||
// I dislike this. But it is necessary if we want to handle those ugly I-P matrices |
|||
inline T convertDistance(state_t tailNode, state_t headNode, T distance) const { |
|||
if (matrixFormat == MatrixFormat::straight) { |
|||
return distance; |
|||
} else { |
|||
if (tailNode == headNode) { |
|||
// diagonal: 1-p = dist |
|||
return one<T>() - distance; |
|||
} else { |
|||
// non-diag: -p = dist |
|||
return zero<T>() - distance; |
|||
} |
|||
} |
|||
} |
|||
|
|||
/*! |
|||
* Returns a map where each state of the input BitVector is mapped to 1 (`one<T>`). |
|||
*/ |
|||
inline StateProbMap allProbOneMap(BitVector bitVector) const { |
|||
StateProbMap stateProbMap; |
|||
for (state_t node : bitVector) { |
|||
stateProbMap.emplace(node, one<T>()); // FIXME check rvalue warning (here and below) |
|||
} |
|||
return stateProbMap; |
|||
} |
|||
|
|||
/*! |
|||
* Given a vector of probabilities so that the `i`th entry corresponds to the |
|||
* probability of state `i`, returns an equivalent map of only the non-zero entries. |
|||
*/ |
|||
inline std::unordered_map<state_t, T> vectorToMap(std::vector<T> probVector) const { |
|||
//assert(probVector.size() == numStates); // numStates may not yet be initialized! // still true? |
|||
|
|||
std::unordered_map<state_t, T> stateProbMap; |
|||
|
|||
for (state_t i = 0; i < probVector.size(); i++) { |
|||
T probEntry = probVector[i]; |
|||
|
|||
// only non-zero entries (i.e. true transitions) are added to the map |
|||
if (probEntry != 0) { |
|||
assert(0 < probEntry <= 1); |
|||
stateProbMap.emplace(i, probEntry); |
|||
} |
|||
} |
|||
|
|||
return stateProbMap; |
|||
} |
|||
|
|||
// ----------------------- |
|||
}; |
|||
} |
|||
} |
|||
} |
|||
|
|||
|
|||
#endif //STORM_UTIL_SHORTESTPATHS_H_ |
|||
@ -0,0 +1,108 @@ |
|||
#include "gtest/gtest.h"
|
|||
#include "storm-config.h"
|
|||
|
|||
#include "storm/builder/ExplicitModelBuilder.h"
|
|||
#include "storm/models/sparse/Dtmc.h"
|
|||
#include "storm/parser/PrismParser.h"
|
|||
#include "storm/storage/SymbolicModelDescription.h"
|
|||
#include "storm/utility/graph.h"
|
|||
#include "storm/utility/shortestPaths.cpp"
|
|||
|
|||
// NOTE: The KSPs / distances 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 ...
|
|||
|
|||
// FIXME: (almost) all of these fail; the question is: is there actually anything wrong or does the new parser yield a different order of states?
|
|||
|
|||
std::shared_ptr<storm::models::sparse::Model<double>> buildExampleModel() { |
|||
std::string prismModelPath = STORM_TEST_RESOURCES_DIR "/dtmc/brp-16-2.pm"; |
|||
storm::storage::SymbolicModelDescription modelDescription = storm::parser::PrismParser::parse(prismModelPath); |
|||
storm::prism::Program program = modelDescription.preprocess().asPrismProgram(); |
|||
return storm::builder::ExplicitModelBuilder<double>(program).build(); |
|||
} |
|||
|
|||
// NOTE: these are hardcoded (obviously), but the model's state indices might change
|
|||
// (e.g., when the parser or model builder are changed)
|
|||
// [state 296 seems to be the new index of the old state 300 (checked a few ksps' probs)]
|
|||
const storm::utility::ksp::state_t testState = 296; |
|||
const storm::utility::ksp::state_t stateWithOnlyOnePath = 1; |
|||
|
|||
TEST(KSPTest, dijkstra) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState); |
|||
|
|||
double dist = spg.getDistance(1); |
|||
EXPECT_DOUBLE_EQ(0.015859334652581887, dist); |
|||
} |
|||
|
|||
TEST(KSPTest, singleTarget) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState); |
|||
|
|||
double dist = spg.getDistance(100); |
|||
EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist); |
|||
} |
|||
|
|||
TEST(KSPTest, reentry) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState); |
|||
|
|||
double dist = spg.getDistance(100); |
|||
EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist); |
|||
|
|||
// get another distance to ensure re-entry is no problem
|
|||
double dist2 = spg.getDistance(500); |
|||
EXPECT_DOUBLE_EQ(3.0462610000679282e-08, dist2); |
|||
} |
|||
|
|||
TEST(KSPTest, groupTarget) { |
|||
auto model = buildExampleModel(); |
|||
auto groupTarget = std::vector<storm::utility::ksp::state_t>{50, 90}; |
|||
auto spg = storm::utility::ksp::ShortestPathsGenerator<double>(*model, groupTarget); |
|||
|
|||
// FIXME comments are outdated (but does it even matter?)
|
|||
// this path should lead to 90
|
|||
double dist1 = spg.getDistance(8); |
|||
EXPECT_DOUBLE_EQ(0.00018449245583999996, dist1); |
|||
|
|||
// this one to 50
|
|||
double dist2 = spg.getDistance(9); |
|||
EXPECT_DOUBLE_EQ(0.00018449245583999996, dist2); |
|||
|
|||
// this one to 90 again
|
|||
double dist3 = spg.getDistance(12); |
|||
EXPECT_DOUBLE_EQ(7.5303043199999984e-06, dist3); |
|||
} |
|||
|
|||
TEST(KSPTest, kTooLargeException) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, stateWithOnlyOnePath); |
|||
|
|||
ASSERT_THROW(spg.getDistance(2), std::invalid_argument); |
|||
} |
|||
|
|||
TEST(KSPTest, kspStateSet) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState); |
|||
|
|||
storm::storage::BitVector referenceBV(model->getNumberOfStates(), false); |
|||
for (auto s : std::vector<storm::utility::ksp::state_t>{0, 1, 3, 5, 7, 10, 14, 19, 25, 29, 33, 36, 40, 44, 50, 56, 62, 70, 77, 85, 92, 98, 104, 112, 119, 127, 134, 140, 146, 154, 161, 169, 176, 182, 188, 196, 203, 211, 218, 224, 230, 238, 245, 253, 260, 266, 272, 281, 288, 296}) { |
|||
referenceBV.set(s, true); |
|||
} |
|||
|
|||
auto bv = spg.getStates(7); |
|||
|
|||
EXPECT_EQ(referenceBV, bv); |
|||
} |
|||
|
|||
TEST(KSPTest, kspPathAsList) { |
|||
auto model = buildExampleModel(); |
|||
storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState); |
|||
|
|||
// TODO: use path that actually has a loop or something to make this more interesting
|
|||
auto reference = storm::utility::ksp::OrderedStateList{296, 288, 281, 272, 266, 260, 253, 245, 238, 230, 224, 218, 211, 203, 196, 188, 182, 176, 169, 161, 154, 146, 140, 134, 127, 119, 112, 104, 98, 92, 85, 77, 70, 62, 56, 50, 44, 36, 29, 40, 33, 25, 19, 14, 10, 7, 5, 3, 1, 0}; |
|||
auto list = spg.getPathAsList(7); |
|||
|
|||
EXPECT_EQ(reference, list); |
|||
} |
|||
Write
Preview
Loading…
Cancel
Save
Reference in new issue