Reimplemented the TopologicalValueIterationNondeterministicLinearEquationSolver with splitting into submatrices.

Added a dtmc example for tests with the StronglyConnectedComponentDecomposition.


Former-commit-id: 0c33793fe6

examples/mdp/scc/scc.pctl

@@ -1,2 +1,2 @@
-Pmin=? [ F end ]
-Pmax=? [ F end ]
+Pmin=? [ (!statetwo) U end ]
+Pmax=? [ (!statetwo) U end ]

src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp

@@ -7,6 +7,11 @@
 #include "src/utility/graph.h"
 #include "src/models/PseudoModel.h"
 #include "src/storage/StronglyConnectedComponentDecomposition.h"
+#include "src/exceptions/IllegalArgumentException.h"
+
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+extern log4cplus::Logger logger;
 
 namespace storm {
     namespace solver {
@@ -39,11 +44,29 @@ namespace storm {
 			//std::vector<std::vector<uint_fast64_t>> stronglyConnectedComponents = storm::utility::graph::performSccDecomposition(this->getModel(), stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
 			//storm::storage::SparseMatrix<T> stronglyConnectedComponentsDependencyGraph = this->getModel().extractSccDependencyGraph(stronglyConnectedComponents);
 
+			std::cout << "TopoSolver Input Matrix: " << A.getRowCount() << " x " << A.getColumnCount() << " with " << A.getEntryCount() << " Entries:" << std::endl;
+
+			uint_fast64_t const rowCount = A.getRowCount();
+			for (uint_fast64_t row = 0; row < rowCount; ++row) {
+				std::cout << "Row " << row << ": ";
+				auto const& rowElement = A.getRow(row);
+				for (auto rowIt = rowElement.begin(); rowIt != rowElement.end(); ++rowIt) {
+					std::cout << rowIt->first << " [" << rowIt->second << "], ";
+				}
+				std::cout << std::endl;
+			}
+
+
 			storm::models::NonDeterministicMatrixBasedPseudoModel<ValueType> pseudoModel(A, nondeterministicChoiceIndices);
-			storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(*static_cast<storm::models::AbstractPseudoModel<ValueType>*>(&pseudoModel), false, false);
-			storm::storage::SparseMatrix<ValueType> stronglyConnectedComponentsDependencyGraph = pseudoModel.extractPartitionDependencyGraph(sccDecomposition);
-			
+			//storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(*static_cast<storm::models::AbstractPseudoModel<ValueType>*>(&pseudoModel), false, false);
+			storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(pseudoModel, false, false);
+
+			if (sccDecomposition.size() == 0) {
+				LOG4CPLUS_ERROR(logger, "Can not solve given Equation System as the SCC Decomposition returned no SCCs.");
+				throw storm::exceptions::IllegalArgumentException() << "Can not solve given Equation System as the SCC Decomposition returned no SCCs.";
+			}
 
+			storm::storage::SparseMatrix<ValueType> stronglyConnectedComponentsDependencyGraph = pseudoModel.extractPartitionDependencyGraph(sccDecomposition);
 			std::vector<uint_fast64_t> topologicalSort = storm::utility::graph::getTopologicalSort(stronglyConnectedComponentsDependencyGraph);
 
 			// Set up the environment for the power method.
@@ -52,12 +75,12 @@ namespace storm {
 				multiplyResult = new std::vector<ValueType>(A.getRowCount());
 				multiplyResultMemoryProvided = false;
 			}
-			std::vector<ValueType>* currentX = &x;
-			bool xMemoryProvided = true;
-			if (newX == nullptr) {
-				newX = new std::vector<ValueType>(x.size());
-				xMemoryProvided = false;
-			}
+			std::vector<ValueType>* currentX = nullptr;
+			//bool xMemoryProvided = true;
+			//if (newX == nullptr) {
+			//	newX = new std::vector<ValueType>(x.size());
+			//	xMemoryProvided = false;
+			//}
 			std::vector<ValueType>* swap = nullptr;
 			uint_fast64_t currentMaxLocalIterations = 0;
 			uint_fast64_t localIterations = 0;
@@ -68,20 +91,62 @@ namespace storm {
 			// solved after all SCCs it depends on have been solved.
 			int counter = 0;
 			std::cout << "Solving Equation System using the TopologicalValueIterationNon..." << std::endl;
+			std::cout << "Found " << sccDecomposition.size() << " SCCs." << std::endl;
+
 			for (auto sccIndexIt = topologicalSort.begin(); sccIndexIt != topologicalSort.end() && converged; ++sccIndexIt) {
 				storm::storage::StateBlock const& scc = sccDecomposition[*sccIndexIt];
 
-				std::cout << "SCC " << counter << " contains:" << std::endl;
-				for (auto sccIt = scc.cbegin(); sccIt != scc.cend(); ++sccIt) {
-					std::cout << *sccIt << ", ";
+				std::cout << "SCC " << counter << " from Index " << *sccIndexIt << " contains:" << std::endl;
+				++counter;
+				for (auto state : scc) {
+					std::cout << state << ", ";
 				}
 				std::cout << std::endl;
 
+				// Generate a submatrix
+				storm::storage::BitVector subMatrixIndices(rowCount, scc.cbegin(), scc.cend());
+				storm::storage::SparseMatrix<ValueType> sccSubmatrix = A.getSubmatrix(subMatrixIndices, nondeterministicChoiceIndices, true);
+				std::vector<ValueType> sccSubB(sccSubmatrix.getRowCount());
+				storm::utility::vector::selectVectorValues<ValueType>(sccSubB, subMatrixIndices, nondeterministicChoiceIndices, b);
+				std::vector<ValueType> sccSubX(sccSubmatrix.getColumnCount());
+				std::vector<ValueType> sccSubXSwap(sccSubmatrix.getColumnCount());
+
+				// Prepare the pointers for swapping in the calculation
+				currentX = &sccSubX;
+				swap = &sccSubXSwap;
+
+				storm::utility::vector::selectVectorValues<ValueType>(sccSubX, subMatrixIndices, x); // x is getCols() large, where as b and multiplyResult are getRows() (nondet. choices times states)
+				std::vector<uint_fast64_t> sccSubNondeterministicChoiceIndices(sccSubmatrix.getColumnCount() + 1);
+				sccSubNondeterministicChoiceIndices.at(0) = 0;
+
+				// Preprocess all dependant states
+				// Remove outgoing transitions and create the ChoiceIndices
+				uint_fast64_t innerIndex = 0;
+				for (uint_fast64_t state: scc) {
+					// Choice Indices
+					sccSubNondeterministicChoiceIndices.at(innerIndex + 1) = sccSubNondeterministicChoiceIndices.at(innerIndex) + (nondeterministicChoiceIndices[state + 1] - nondeterministicChoiceIndices[state]);
+
+					for (auto rowGroupIt = nondeterministicChoiceIndices[state]; rowGroupIt != nondeterministicChoiceIndices[state + 1]; ++rowGroupIt) {
+						storm::storage::SparseMatrix<ValueType>::const_rows row = A.getRow(rowGroupIt);
+						for (auto rowIt = row.begin(); rowIt != row.end(); ++rowIt) {
+							if (!subMatrixIndices.get(rowIt->first)) {
+								// This is an outgoing transition of a state in the SCC to a state not included in the SCC
+								// Subtracting Pr(tau) * x_other from b fixes that
+								sccSubB.at(innerIndex) = sccSubB.at(innerIndex) - (rowIt->second * x.at(rowIt->first));
+							}
+						}
+					}
+					++innerIndex;
+				}
+
 				// For the current SCC, we need to perform value iteration until convergence.
 				localIterations = 0;
 				converged = false;
 				while (!converged && localIterations < this->maximalNumberOfIterations) {
 					// Compute x' = A*x + b.
+					sccSubmatrix.multiplyWithVector(*currentX, *multiplyResult);
+					storm::utility::vector::addVectorsInPlace<ValueType>(*multiplyResult, sccSubB);
+
 					//A.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
 					//storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
 
@@ -93,43 +158,46 @@ namespace storm {
 
 					// Reduce the vector x' by applying min/max for all non-deterministic choices.
 					if (minimize) {
-						//storm::utility::reduceVectorMin(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
+						storm::utility::vector::reduceVectorMin<ValueType>(*multiplyResult, *swap, sccSubNondeterministicChoiceIndices);
 					}
 					else {
-						//storm::utility::reduceVectorMax(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
+						storm::utility::vector::reduceVectorMax<ValueType>(*multiplyResult, *swap, sccSubNondeterministicChoiceIndices);
 					}
 
 					// Determine whether the method converged.
 					// TODO: It seems that the equalModuloPrecision call that compares all values should have a higher
 					// running time. In fact, it is faster. This has to be investigated.
 					// converged = storm::utility::equalModuloPrecision(*currentX, *newX, scc, precision, relative);
-					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, this->precision, this->relative);
+					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *swap, this->precision, this->relative);
 
 					// Update environment variables.
-					swap = currentX;
-					currentX = newX;
-					newX = swap;
+					std::swap(currentX, swap);
+
 					++localIterations;
 					++globalIterations;
 				}
 
+				// The Result of this SCC has to be taken back into the main result vector
+				innerIndex = 0;
+				for (uint_fast64_t state: scc) {
+					x.at(state) = currentX->at(innerIndex);
+					++innerIndex;
+				}
+
+				// Since the pointers for swapping in the calculation point to temps they should not be valide anymore
+				currentX = nullptr;
+				swap = nullptr;
+
 				// As the "number of iterations" of the full method is the maximum of the local iterations, we need to keep
 				// track of the maximum.
 				if (localIterations > currentMaxLocalIterations) {
 					currentMaxLocalIterations = localIterations;
 				}
 			}
-
-			// If we performed an odd number of global iterations, we need to swap the x and currentX, because the newest
-			// result is currently stored in currentX, but x is the output vector.
-			// TODO: Check whether this is correct or should be put into the for-loop over SCCs.
-			if (globalIterations % 2 == 1) {
-				std::swap(x, *currentX);
-			}
 			
-			if (!xMemoryProvided) {
-				delete newX;
-			}
+			//if (!xMemoryProvided) {
+			//	delete newX;
+			//}
 
 			if (!multiplyResultMemoryProvided) {
 				delete multiplyResult;

src/storage/StronglyConnectedComponentDecomposition.cpp

@@ -35,7 +35,7 @@ namespace storm {
         template <typename ValueType>
         StronglyConnectedComponentDecomposition<ValueType>::StronglyConnectedComponentDecomposition(storm::models::AbstractModel<ValueType> const& model, storm::storage::BitVector const& subsystem, bool dropNaiveSccs, bool onlyBottomSccs) {
 			storm::models::ModelBasedPseudoModel<ValueType> encapsulatedModel(model);
-			performSccDecomposition(*static_cast<storm::models::AbstractPseudoModel<ValueType>*>(&encapsulatedModel), subsystem, dropNaiveSccs, onlyBottomSccs);
+			performSccDecomposition(encapsulatedModel, subsystem, dropNaiveSccs, onlyBottomSccs);
         }
 
 		template <typename ValueType>

src/storm.cpp

@@ -17,6 +17,7 @@
 #include <iostream>
 #include <fstream>
 #include <cstdio>
+#include <climits>
 #include <sstream>
 #include <vector>
 
@@ -29,6 +30,7 @@
 #include "src/models/AtomicPropositionsLabeling.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
+#include "src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
 #include "src/solver/GurobiLpSolver.h"
@@ -132,6 +134,16 @@ void setUpFileLogging() {
 	logger.addAppender(fileLogAppender);
 }
 
+/*!
+* Gives the current working directory
+*
+* @return std::string The path of the current working directory
+*/
+std::string getCurrentWorkingDirectory() {
+	char temp[512];
+	return (getcwd(temp, 512 - 1) ? std::string(temp) : std::string(""));
+}
+
 /*!
  * Prints the header.
  */
@@ -146,7 +158,8 @@ void printHeader(const int argc, const char* argv[]) {
 	for (int i = 0; i < argc; ++i) {
 		commandStream << argv[i] << " ";
 	}
-	std::cout << "Command line: " << commandStream.str() << std::endl << std::endl;
+	std::cout << "Command line: " << commandStream.str() << std::endl;
+	std::cout << "Current working directory: " << getCurrentWorkingDirectory() << std::endl << std::endl;
 }
 
 /*!
@@ -234,7 +247,8 @@ storm::modelchecker::prctl::AbstractModelChecker<double>* createPrctlModelChecke
  */
 storm::modelchecker::prctl::AbstractModelChecker<double>* createPrctlModelChecker(storm::models::Mdp<double>& mdp) {
     // Create the appropriate model checker.
-    return new storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double>(mdp);
+    //return new storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double>(mdp);
+	return new storm::modelchecker::prctl::TopologicalValueIterationMdpPrctlModelChecker<double>(mdp);
 }
 
 /*!

src/utility/OsDetection.h

@@ -4,10 +4,12 @@
 #if defined __linux__ || defined __linux
 #	define LINUX
 #	include <sys/mman.h>
+#	include <unistd.h>
 #include <execinfo.h> // Required by ErrorHandling.h
 #include <cxxabi.h> // Required by ErrorHandling.h
 #include <sys/time.h> // Required by storm.cpp, Memory Usage
 #include <sys/resource.h> // Required by storm.cpp, Memory Usage
+#	define GetCurrentDir getcwd
 #elif defined TARGET_OS_MAC || defined __apple__ || defined __APPLE__
 #	define MACOSX
 #	define _DARWIN_USE_64_BIT_INODE
@@ -17,6 +19,7 @@
 #	include <cxxabi.h> // Required by ErrorHandling.h
 #	include <sys/time.h> // Required by storm.cpp, Memory Usage
 #	include <sys/resource.h> // Required by storm.cpp, Memory Usage
+#	define GetCurrentDir getcwd
 #elif defined _WIN32 || defined _WIN64
 #	define WINDOWS
 #	ifndef NOMINMAX
@@ -28,8 +31,10 @@
 #	include <winnt.h>
 #	include <DbgHelp.h>
 #	include <Psapi.h>
+#	include <direct.h>
 #	define strncpy strncpy_s
 #	define sscanf sscanf_s
+#	define GetCurrentDir _getcwd
 
 // This disables Warning C4250 - Diamond Inheritance Dominance
 #pragma warning(disable:4250)

test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp

@@ -7,6 +7,36 @@
 #include "src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h"
 #include "src/parser/AutoParser.h"
 
+TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, SmallLinEqSystem) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder(3, 3);
+	ASSERT_NO_THROW(matrixBuilder.addNextValue(0, 2, 1.0));
+	ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 0, 4.0));
+	ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 1, 7.0));
+	ASSERT_NO_THROW(matrixBuilder.addNextValue(2, 1, 7.0));
+	ASSERT_NO_THROW(matrixBuilder.addNextValue(2, 2, -1.0));
+
+	storm::storage::SparseMatrix<double> matrix;
+	ASSERT_NO_THROW(matrix = matrixBuilder.build());
+
+	ASSERT_EQ(3, matrix.getRowCount());
+	ASSERT_EQ(3, matrix.getColumnCount());
+	ASSERT_EQ(5, matrix.getEntryCount());
+	
+	// Solve the Linear Equation System
+	storm::solver::TopologicalValueIterationNondeterministicLinearEquationSolver<double> topoSolver;
+
+	std::vector<double> x(3);
+	std::vector<double> b = { 5, 8, 2 };
+	std::vector<uint_fast64_t> choices = { 0, 1, 2, 3 };
+
+	ASSERT_NO_THROW(topoSolver.solveEquationSystem(true, matrix, x, b, choices));
+
+	storm::settings::Settings* s = storm::settings::Settings::getInstance();
+	ASSERT_LT(std::abs(x.at(0) - 0.25), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+	ASSERT_LT(std::abs(x.at(1) - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+	ASSERT_LT(std::abs(x.at(2) - 5.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
 TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	storm::settings::Settings* s = storm::settings::Settings::getInstance();
 	//storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
@@ -17,7 +47,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
     
 	ASSERT_EQ(mdp->getNumberOfStates(), 11ull);
-	ASSERT_EQ(mdp->getNumberOfTransitions(), 17ull);
+	ASSERT_EQ(mdp->getNumberOfTransitions(), 18ull);
     
 	storm::modelchecker::prctl::TopologicalValueIterationMdpPrctlModelChecker<double> mc(*mdp);
     

test/functional/storage/StronglyConnectedComponentDecompositionTest.cpp

@@ -44,3 +44,31 @@ TEST(StronglyConnectedComponentDecomposition, FullSystem2) {
 
     markovAutomaton = nullptr;
 }
+
+TEST(StronglyConnectedComponentDecomposition, MatrixBasedSystem) {
+	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/dtmc/scc/scc.tra", STORM_CPP_BASE_PATH "/examples/dtmc/scc/scc.lab", "", "");
+	std::shared_ptr<storm::models::Dtmc<double>> dtmc = parser.getModel<storm::models::Dtmc<double>>();
+
+	storm::storage::StronglyConnectedComponentDecomposition<double> sccDecomposition;
+	ASSERT_NO_THROW(sccDecomposition = storm::storage::StronglyConnectedComponentDecomposition<double>(*dtmc, true, false));
+
+	ASSERT_EQ(sccDecomposition.size(), 3);
+
+	// Now, because there is no ordering we have to check the contents of the MECs in a symmetrical way.
+	storm::storage::StateBlock const& scc1 = sccDecomposition[0];
+	storm::storage::StateBlock const& scc2 = sccDecomposition[1];
+	storm::storage::StateBlock const& scc3 = sccDecomposition[2];
+
+	std::vector<uint_fast64_t> correctScc1 = { 1, 2, 3, 4 };
+	std::vector<uint_fast64_t> correctScc2 = { 5, 6, 7, 8 };
+	std::vector<uint_fast64_t> correctScc3 = { 0 };
+
+	ASSERT_TRUE(scc1 == storm::storage::StateBlock(correctScc1.begin(), correctScc1.end()) || scc1 == storm::storage::StateBlock(correctScc2.begin(), correctScc2.end()) || scc1 == storm::storage::StateBlock(correctScc3.begin(), correctScc3.end()));
+	ASSERT_TRUE(scc2 == storm::storage::StateBlock(correctScc1.begin(), correctScc1.end()) || scc2 == storm::storage::StateBlock(correctScc2.begin(), correctScc2.end()) || scc2 == storm::storage::StateBlock(correctScc3.begin(), correctScc3.end()));
+	ASSERT_TRUE(scc3 == storm::storage::StateBlock(correctScc1.begin(), correctScc1.end()) || scc3 == storm::storage::StateBlock(correctScc2.begin(), correctScc2.end()) || scc3 == storm::storage::StateBlock(correctScc3.begin(), correctScc3.end()));
+
+	ASSERT_NO_THROW(sccDecomposition = storm::storage::StronglyConnectedComponentDecomposition<double>(*dtmc, true, true));
+	ASSERT_EQ(2, sccDecomposition.size());
+
+	dtmc = nullptr;
+}