From fe6804e16426e80b9a718f091a24210e7d5ba835 Mon Sep 17 00:00:00 2001
From: tomjanson <tom.janson@rwth-aachen.de>
Date: Fri, 18 Nov 2016 16:52:07 +0100
Subject: [PATCH] KSP: matrix format conversion & lots of type stuff

# Conflicts:
#	src/python/storm-tom.cpp
#	src/utility/shortestPaths.cpp
#	src/utility/shortestPaths.h
#	test/functional/utility/KSPTest.cpp
#	test/functional/utility/PdtmcInstantiationTest.cpp
---
 src/utility/shortestPaths.cpp       | 32 +++++++------
 src/utility/shortestPaths.h         | 69 ++++++++++++++++++++++-------
 test/functional/utility/KSPTest.cpp | 16 +++----
 3 files changed, 79 insertions(+), 38 deletions(-)

diff --git a/src/utility/shortestPaths.cpp b/src/utility/shortestPaths.cpp
index dda452d5a..ad25e2616 100644
--- a/src/utility/shortestPaths.cpp
+++ b/src/utility/shortestPaths.cpp
@@ -7,12 +7,16 @@ 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) :
+            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) {
+                    targetProbMap(targetProbMap),
+                    matrixFormat(matrixFormat) {
 
                 computePredecessors();
 
@@ -28,28 +32,29 @@ namespace storm {
             }
 
             template <typename T>
-            ShortestPathsGenerator<T>::ShortestPathsGenerator(storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& targetProbVector, BitVector const& initialStates)
-                    : ShortestPathsGenerator<T>(transitionMatrix, vectorToMap(targetProbVector), initialStates) {}
+            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(std::shared_ptr<models::sparse::Model<T>> model, BitVector const& targetBV)
-                    : ShortestPathsGenerator<T>(model->getTransitionMatrix(), allProbOneMap(targetBV), model->getInitialStates()) {}
+            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(std::shared_ptr<models::sparse::Model<T>> model, state_t singleTarget)
+            ShortestPathsGenerator<T>::ShortestPathsGenerator(Model const& 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, std::vector<state_t> const& targetList)
-                    : ShortestPathsGenerator<T>(model, BitVector(model->getNumberOfStates(), targetList)) {}
+            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(std::shared_ptr<models::sparse::Model<T>> model, std::string const& targetLabel)
-                    : ShortestPathsGenerator<T>(model, model->getStates(targetLabel)) {}
+            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) {
@@ -152,7 +157,8 @@ namespace storm {
                             state_t otherNode = transition.getColumn();
 
                             // note that distances are probabilities, thus they are multiplied and larger is better
-                            T alternateDistance = shortestPathDistances[currentNode] * transition.getValue();
+                            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);
@@ -225,7 +231,7 @@ namespace storm {
 
                     for (auto const& transition : transitionMatrix.getRowGroup(tailNode)) {
                         if (transition.getColumn() == headNode) {
-                            return transition.getValue();
+                            return convertDistance(tailNode, headNode, transition.getValue());
                         }
                     }
 
diff --git a/src/utility/shortestPaths.h b/src/utility/shortestPaths.h
index 67e857e5d..fb1e0303c 100644
--- a/src/utility/shortestPaths.h
+++ b/src/utility/shortestPaths.h
@@ -23,9 +23,11 @@
 namespace storm {
     namespace utility {
         namespace ksp {
-            typedef storage::sparse::state_type state_t;
+            using state_t = storage::sparse::state_type;
+            using OrderedStateList = std::vector<state_t>;
             using BitVector = storage::BitVector;
-            typedef std::vector<state_t> ordered_state_list_t;
+
+            // -- helper structs/classes -----------------------------------------------------------------------------
 
             template <typename T>
             struct Path {
@@ -46,29 +48,38 @@ namespace storm {
                 }
             };
 
+            // 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(std::shared_ptr<models::sparse::Model<T>> model, BitVector const& targetBV);
+                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(std::shared_ptr<models::sparse::Model<T>> model, state_t singleTarget);
-                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, std::vector<state_t> const& targetList);
-                ShortestPathsGenerator(std::shared_ptr<models::sparse::Model<T>> model, std::string const& targetLabel = "target");
+                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(storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& targetProbVector, BitVector const& initialStates);
-                ShortestPathsGenerator(storage::SparseMatrix<T> const& maybeTransitionMatrix, std::unordered_map<state_t, T> const& targetProbMap, BitVector const& initialStates);
+                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(){}
@@ -93,19 +104,21 @@ namespace storm {
                  * Computes KSP if not yet computed.
                  * @throws std::invalid_argument if no such k-shortest path exists
                  */
-                ordered_state_list_t getPathAsList(unsigned long k);
+                OrderedStateList getPathAsList(unsigned long k);
 
 
             private:
-                storage::SparseMatrix<T> const& transitionMatrix; // FIXME: store reference instead (?)
+                Matrix const& transitionMatrix; // FIXME: store reference instead (?)
                 state_t numStates; // includes meta-target, i.e. states in model + 1
                 state_t metaTarget;
                 BitVector initialStates;
-                std::unordered_map<state_t, T> targetProbMap;
+                StateProbMap targetProbMap;
+
+                MatrixFormat matrixFormat;
 
-                std::vector<ordered_state_list_t>     graphPredecessors;
+                std::vector<OrderedStateList>         graphPredecessors;
                 std::vector<boost::optional<state_t>> shortestPathPredecessors;
-                std::vector<ordered_state_list_t>     shortestPathSuccessors;
+                std::vector<OrderedStateList>         shortestPathSuccessors;
                 std::vector<T>                        shortestPathDistances;
 
                 std::vector<std::vector<Path<T>>> kShortestPaths;
@@ -163,6 +176,7 @@ namespace storm {
 
 
                 // --- tiny helper fcts ---
+
                 inline bool isInitialState(state_t node) const {
                     return find(initialStates.begin(), initialStates.end(), node) != initialStates.end();
                 }
@@ -171,25 +185,45 @@ namespace storm {
                     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 std::unordered_map<state_t, T> allProbOneMap(BitVector bitVector) const {
-                    std::unordered_map<state_t, T> stateProbMap;
+                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);
 
                     std::unordered_map<state_t, T> stateProbMap;
 
-                    // non-zero entries (i.e. true transitions) are added to the map
                     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);
@@ -198,6 +232,7 @@ namespace storm {
 
                     return stateProbMap;
                 }
+
                 // -----------------------
             };
         }
diff --git a/test/functional/utility/KSPTest.cpp b/test/functional/utility/KSPTest.cpp
index d1c9a2f8f..2389f3f1c 100644
--- a/test/functional/utility/KSPTest.cpp
+++ b/test/functional/utility/KSPTest.cpp
@@ -26,7 +26,7 @@ TEST(KSPTest, dijkstra) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 300;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState);
 
     double dist = spg.getDistance(1);
     EXPECT_DOUBLE_EQ(0.015859334652581887, dist);
@@ -36,7 +36,7 @@ TEST(KSPTest, singleTarget) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 300;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState);
 
     double dist = spg.getDistance(100);
     EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist);
@@ -46,7 +46,7 @@ TEST(KSPTest, reentry) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 300;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState);
 
     double dist = spg.getDistance(100);
     EXPECT_DOUBLE_EQ(1.5231305000339649e-06, dist);
@@ -60,7 +60,7 @@ 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);
+    auto spg = storm::utility::ksp::ShortestPathsGenerator<double>(*model, groupTarget);
 
     // this path should lead to 90
     double dist1 = spg.getDistance(8);
@@ -79,7 +79,7 @@ TEST(KSPTest, kTooLargeException) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 1;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    storm::utility::ksp::ShortestPathsGenerator<double> spg(*model, testState);
 
     ASSERT_THROW(spg.getDistance(2), std::invalid_argument);
 }
@@ -88,7 +88,7 @@ TEST(KSPTest, kspStateSet) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 300;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    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>{300, 293, 285, 279, 271, 265, 259, 252, 244, 237, 229, 223, 217, 210, 202, 195, 187, 181, 175, 168, 160, 153, 145, 139, 133, 126, 118, 111, 103, 97, 91, 84, 76, 69, 61, 55, 49, 43, 35, 41, 32, 25, 19, 14, 10, 7, 4, 2, 1, 0}) {
@@ -104,10 +104,10 @@ TEST(KSPTest, kspPathAsList) {
 	auto model = buildExampleModel();
 
     storm::utility::ksp::state_t testState = 300;
-    storm::utility::ksp::ShortestPathsGenerator<double> spg(model, testState);
+    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::ordered_state_list_t{300, 293, 285, 279, 271, 265, 259, 252, 244, 237, 229, 223, 217, 210, 202, 195, 187, 181, 175, 168, 160, 153, 145, 139, 133, 126, 118, 111, 103, 97, 91, 84, 76, 69, 61, 55, 49, 43, 35, 41, 32, 25, 19, 14, 10, 7, 4, 2, 1, 0};
+    auto reference = storm::utility::ksp::OrderedStateList{300, 293, 285, 279, 271, 265, 259, 252, 244, 237, 229, 223, 217, 210, 202, 195, 187, 181, 175, 168, 160, 153, 145, 139, 133, 126, 118, 111, 103, 97, 91, 84, 76, 69, 61, 55, 49, 43, 35, 41, 32, 25, 19, 14, 10, 7, 4, 2, 1, 0};
     auto list = spg.getPathAsList(7);
 
     EXPECT_EQ(list, reference);