Merge branch 'critSubsysMerge'

Conflicts:
	src/utility/StormOptions.cpp

Former-commit-id: 924ecb35f7

src/counterexamples/PathBasedSubsystemGenerator.h

@@ -0,0 +1,716 @@
+/*
+ * PathBasedSubsystemGenerator.h
+ *
+ *  Created on: 11.10.2013
+ *      Author: Manuel Sascha Weiand
+ */
+
+#ifndef STORM_COUNTEREXAMPLES_PATHBASEDSUBSYSTEMGENERATOR_H_
+#define STORM_COUNTEREXAMPLES_PATHBASEDSUBSYSTEMGENERATOR_H_
+
+#include "src/models/Dtmc.h"
+#include "src/models/AbstractModel.h"
+#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "src/solver/GmmxxLinearEquationSolver.h"
+#include "src/storage/BitVector.h"
+#include "src/storage/SparseMatrix.h"
+
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
+extern log4cplus::Logger logger;
+
+namespace storm {
+    namespace counterexamples {
+
+    template <class T>
+    class PathBasedSubsystemGenerator {
+
+private:
+
+	template <class CompareType>
+	class CompareStates {
+	public:
+		bool operator()(const std::pair<uint_fast64_t, CompareType>& s1, const std::pair<uint_fast64_t, CompareType>& s2) {
+			return s1.second < s2.second;
+		}
+	};
+
+public:
+
+	/*!
+	 *
+	 */
+	static void computeShortestDistances(storm::storage::SparseMatrix<T> const& transMat, storm::storage::BitVector& subSysStates, storm::storage::BitVector& terminalStates, storm::storage::BitVector& allowedStates, std::vector<std::pair<uint_fast64_t, T>>& distances) {
+
+		std::multiset<std::pair<uint_fast64_t, T>, CompareStates<T> > activeSet;
+		//std::priority_queue<std::pair<uint_fast64_t, T*>, std::vector<std::pair<uint_fast64_t, T*>>, CompareStates<T> > queue;
+
+		// resize and init distances
+		const std::pair<uint_fast64_t, T> initDistances(0, (T) -1);
+		distances.resize(transMat.getColumnCount(), initDistances);
+
+		//since gcc 4.7.2 does not implement std::set::emplace(), insert is used.
+		std::pair<uint_fast64_t, T> state;
+
+		// First store all transitions from initial states
+		// Also save all found initial states in array of discovered states.
+		for(storm::storage::BitVector::constIndexIterator init = subSysStates.begin(); init != subSysStates.end(); ++init) {
+			//use init state only if it is allowed
+			if(allowedStates.get(*init)) {
+
+				if(terminalStates.get(*init)) {
+						// it's an init -> target search
+						// save target state as discovered and get it's outgoing transitions
+
+						distances[*init].first = *init;
+						distances[*init].second = (T) 1;
+				}
+
+				typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIt = transMat.begin(*init);
+				for(auto trans = rowIt.begin() ; trans != rowIt.end(); ++trans) {
+					//save transition only if it's no 'virtual transition of prob 0 and it doesn't go from init state to init state.
+					if(trans.value() != (T) 0 && !subSysStates.get(trans.column())) {
+						//new state?
+						if(distances[trans.column()].second == (T) -1) {
+							distances[trans.column()].first = *init;
+							distances[trans.column()].second = trans.value();
+
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+						}
+						else if(distances[trans.column()].second < trans.value()){
+							//This state has already been discovered
+							//And the distance can be improved by using this transition.
+
+							//find state in set, remove it, reenter it with new and correct values.
+							auto range = activeSet.equal_range(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+							for(;range.first != range.second; range.first++) {
+								if(trans.column() == range.first->first) {
+									activeSet.erase(range.first);
+									break;
+								}
+							}
+
+							distances[trans.column()].first = *init;
+							distances[trans.column()].second = trans.value();
+
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), trans.value()));
+						}
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Initialized.");
+
+		//Now find the shortest distances to all states
+		while(!activeSet.empty()) {
+
+			// copy here since using a reference leads to segfault
+			std::pair<uint_fast64_t, T> activeState = *(--activeSet.end());
+			activeSet.erase(--activeSet.end());
+
+			// If this is an initial state, do not consider its outgoing transitions, since all relevant ones have already been considered
+			// Same goes for forbidden states since they may not be used on a path, except as last node.
+			if(!subSysStates.get(activeState.first) && allowedStates.get(activeState.first)) {
+				// Look at all neighbors
+				typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIt = transMat.begin(activeState.first);
+				for(typename storm::storage::SparseMatrix<T>::ConstIterator trans = rowIt.begin(); trans != rowIt.end(); ++trans) {
+					// Only consider the transition if it's not virtual
+					if(trans.value() != (T) 0) {
+
+						T distance = activeState.second * trans.value();
+
+						//not discovered or initial terminal state
+						if(distances[trans.column()].second == (T)-1) {
+							//New state discovered -> save it
+							distances[trans.column()].first = activeState.first;
+							distances[trans.column()].second = distance;
+
+							// push newly discovered state into activeSet
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
+						}
+						else if(distances[trans.column()].second < distance ){
+							//This state has already been discovered
+							//And the distance can be improved by using this transition.
+
+							//find state in set, remove it, reenter it with new and correct values.
+
+							auto range = activeSet.equal_range(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+							for(;range.first != range.second; range.first++) {
+								if(trans.column() == range.first->first) {
+									activeSet.erase(range.first);
+									break;
+								}
+							}
+
+							distances[trans.column()].first = activeState.first;
+							distances[trans.column()].second = distance;
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
+						}
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Discovery done.");
+	}
+
+	/*!
+	 *
+	 */
+	static void doDijkstraSearch(storm::storage::SparseMatrix<T> const& transMat, storm::storage::BitVector& subSysStates, storm::storage::BitVector& terminalStates, storm::storage::BitVector& allowedStates, std::vector<std::pair<uint_fast64_t, T>>& itDistances, std::vector<std::pair<uint_fast64_t, T>>& distances) {
+		std::multiset<std::pair<uint_fast64_t, T>, CompareStates<T> > activeSet;
+
+		// resize and init distances
+		const std::pair<uint_fast64_t, T> initDistances(0, (T) -1);
+		distances.resize(transMat.getColumnCount(), initDistances);
+
+		//since gcc 4.7.2 does not implement std::set::emplace(), insert is used.
+		std::pair<uint_fast64_t, T> state;
+
+		// First store all transitions from initial states
+		// Also save all found initial states in array of discovered states.
+		for(storm::storage::BitVector::constIndexIterator init = subSysStates.begin(); init != subSysStates.end(); ++init) {
+			//use init state only if it is allowed
+			if(allowedStates.get(*init)) {
+
+				if(terminalStates.get(*init)) {
+						// it's a subsys -> subsys search
+						// ignore terminal state completely
+						// (since any target state that is only reached by a path going through this state should not be reached)
+						continue;
+				}
+
+				typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIt = transMat.begin(*init);
+				for(typename storm::storage::SparseMatrix<T>::ConstIterator trans = rowIt.begin(); trans != rowIt.end(); ++trans) {
+					//save transition only if it's no 'virtual transition of prob 0 and it doesn't go from init state to init state.
+					if(trans.value() != (T) 0 && !subSysStates.get(trans.column())) {
+						//new state?
+						if(distances[trans.column()].second == (T) -1) {
+							//for initialization of subsys -> subsys search use prob (init -> subsys state -> found state) instead of prob(subsys state -> found state)
+							distances[trans.column()].first = *init;
+							distances[trans.column()].second = trans.value() * (itDistances[*init].second == -1 ? 1 : itDistances[*init].second);
+
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+						}
+						else if(distances[trans.column()].second < trans.value() * itDistances[*init].second){
+							//This state has already been discovered
+							//And the distance can be improved by using this transition.
+
+							//find state in set, remove it, reenter it with new and correct values.
+							auto range = activeSet.equal_range(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+							for(;range.first != range.second; range.first++) {
+								if(trans.column() == range.first->first) {
+									activeSet.erase(range.first);
+									break;
+								}
+							}
+
+							//for initialization of subsys -> subsys search use prob (init -> subsys state -> found state) instead of prob(subsys state -> found state)
+							distances[trans.column()].first = *init;
+							distances[trans.column()].second = trans.value() * (itDistances[*init].second == -1 ? 1 : itDistances[*init].second);
+
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), trans.value()));
+						}
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Initialized.");
+
+		//Now find the shortest distances to all states
+		while(!activeSet.empty()) {
+
+			// copy here since using a reference leads to segfault
+			std::pair<uint_fast64_t, T> activeState = *(--activeSet.end());
+			activeSet.erase(--activeSet.end());
+
+			// Always stop at first target/terminal state
+			//if(terminalStates.get(activeState.first) || subSysStates.get(activeState.first)) break;
+
+			// If this is an initial state, do not consider its outgoing transitions, since all relevant ones have already been considered
+			// Same goes for forbidden states since they may not be used on a path, except as last node.
+			if(!subSysStates.get(activeState.first) && allowedStates.get(activeState.first)) {
+				// Look at all neighbors
+				typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIt = transMat.begin(activeState.first);
+				for(typename storm::storage::SparseMatrix<T>::ConstIterator trans = rowIt.begin(); trans != rowIt.end(); ++trans) {
+					// Only consider the transition if it's not virtual
+					if(trans.value() != (T) 0) {
+
+						T distance = activeState.second * trans.value();
+
+						//not discovered or initial terminal state
+						if(distances[trans.column()].second == (T)-1) {
+							//New state discovered -> save it
+							distances[trans.column()].first = activeState.first;
+							distances[trans.column()].second = distance;
+
+							// push newly discovered state into activeSet
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
+						}
+						else if(distances[trans.column()].second < distance ){
+							//This state has already been discovered
+							//And the distance can be improved by using this transition.
+
+							//find state in set, remove it, reenter it with new and correct values.
+
+							auto range = activeSet.equal_range(std::pair<uint_fast64_t, T>(trans.column(), distances[trans.column()].second));
+							for(;range.first != range.second; range.first++) {
+								if(trans.column() == range.first->first) {
+									activeSet.erase(range.first);
+									break;
+								}
+							}
+
+							distances[trans.column()].first = activeState.first;
+							distances[trans.column()].second = distance;
+							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
+						}
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Discovery done.");
+	}
+
+	/*!
+	 *
+	 *//*
+	template <typename T>
+	static void doBackwardsSearch(storm::storage::SparseMatrix<T> const& transMat, storm::storage::BitVector& initStates, storm::storage::BitVector& subSysStates, storm::storage::BitVector& terminalStates, storm::storage::BitVector& allowedStates, std::vector<T>& probabilities, std::vector<uint_fast64_t>& shortestPath, T& probability) {
+		std::multiset<std::pair<uint_fast64_t, T>, CompareStates<T> > activeSet;
+
+		// resize and init distances
+		const std::pair<uint_fast64_t, T> initDistances(0, (T) -1);
+		std::vector<std::pair<uint_fast64_t, T>> distances(transMat.getColumnCount(), initDistances);
+
+		//since the transition matrix only gives a means to iterate over successors and not over predecessors and there is no Transpose for the matrix
+		//GraphTransitions is used
+		storm::models::GraphTransitions<T> backTrans(transMat, false);
+
+		//First store all allowed predecessors of target states that are not in the subsystem
+		for(storm::storage::BitVector::constIndexIterator target = terminalStates.begin(); target != terminalStates.end(); ++target) {
+
+			// if there is a terminal state that is an initial state then prob == 1 and return
+			if(initStates.get(*target)){
+				distances[*target].first = *target;
+				distances[*target].second = (T) 1;
+				return;
+			}
+
+			//iterate over predecessors
+			for(auto iter = backTrans.beginStateSuccessorsIterator(*target); iter != backTrans.endStateSuccessorsIterator(*target); iter++) {
+				//only use if allowed and not in subsys and not terminal
+				if(allowedStates.get(*iter) && !subSysStates.get(*iter) && !terminalStates.get(*iter)) {
+					//new state?
+					if(distances[*iter].second == (T) -1) {
+						// save as discovered and push into active set
+						distances[*iter].first = *target;												//successor
+						distances[*iter].second = transMat.getValue(*iter, *target);					//prob of shortest path
+
+						activeSet.insert(std::pair<uint_fast64_t, T>(*iter, probabilities[*iter]));		//prob of reaching some terminal state from pred.
+					}
+					else {
+						// state was already discovered
+						// is this the better transition?
+						if(distances[*iter].second > transMat.getValue(*iter, *target)) {
+							distances[*iter].first = *target;
+							distances[*iter].second = transMat.getValue(*iter, *target);
+						}
+					}
+				}
+			}
+		}
+
+		//Now store all allowed predecessors of subsystem states that are not subsystem states themselves
+		for(storm::storage::BitVector::constIndexIterator sysState = subSysStates.begin(); sysState != subSysStates.end(); ++sysState) {
+			//iterate over predecessors
+			for(auto iter = backTrans.beginStateSuccessorsIterator(*sysState); iter != backTrans.endStateSuccessorsIterator(*sysState); iter++) {
+				//only use if allowed and not in subsys and not terminal
+				if(allowedStates.get(*iter) && !subSysStates.get(*iter) && !terminalStates.get(*iter)) {
+					//new state?
+					if(distances[*iter].second == (T) -1) {
+						// save as discovered and push into active set
+						distances[*iter].first = *sysState;												//successor
+						distances[*iter].second = transMat.getValue(*iter, *sysState);					//prob of shortest path
+
+						activeSet.insert(std::pair<uint_fast64_t, T>(*iter, probabilities[*iter]));		//prob of reaching some terminal state from pred.
+					}
+					else {
+						// state was already discovered
+						// is this the better transition?
+						if(distances[*iter].second > transMat.getValue(*iter, *sysState)) {
+							distances[*iter].first = *sysState;
+							distances[*iter].second = transMat.getValue(*iter, *sysState);
+						}
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Initialized.");
+
+		// Do the backwards search
+		std::pair<uint_fast64_t, T> state;
+		uint_fast64_t activeState;
+		while(!activeSet.empty()) {
+			// copy here since using a reference leads to segfault
+			state = *(--activeSet.end());
+			activeState = state.first;
+			activeSet.erase(--activeSet.end());
+
+			//stop on the first subsys/init state
+			if(initStates.get(activeState) || subSysStates.get(activeState)) break;
+
+			// If this is a subSys or terminal state, do not consider its incoming transitions, since all relevant ones have already been considered
+			if(!terminalStates.get(activeState) && !subSysStates.get(activeState)) {
+				//iterate over predecessors
+				for(auto iter = backTrans.beginStateSuccessorsIterator(activeState); iter != backTrans.endStateSuccessorsIterator(activeState); iter++) {
+					//only if transition is not "virtual" and no selfloop
+					if(*iter != activeState && transMat.getValue(*iter, activeState) != (T) 0) {
+						//new state?
+						if(distances[*iter].second == (T) -1) {
+							// save as discovered and push into active set
+							distances[*iter].first = activeState;
+							distances[*iter].second = transMat.getValue(*iter, activeState) * distances[activeState].second;
+
+							activeSet.insert(std::pair<uint_fast64_t, T>(*iter, probabilities[*iter]));
+						}
+						else {
+							// state was already discovered
+							// is this the better transition?
+							if(distances[*iter].second < transMat.getValue(*iter, activeState) * distances[activeState].second) {
+								distances[*iter].first = activeState;
+								distances[*iter].second = transMat.getValue(*iter, activeState) * distances[activeState].second;
+							}
+						}
+					}
+				}
+			}
+		}
+
+		//get path probability
+		probability = distances[activeState].second;
+		if(probability == (T) -1) probability = 1;
+
+		// iterate over the successors until reaching the end of the finite path
+		shortestPath.push_back(activeState);
+		activeState = distances[activeState].first;
+		while(!terminalStates.get(activeState) && !subSysStates.get(activeState)) {
+			shortestPath.push_back(activeState);
+			activeState = distances[activeState].first;
+		}
+		shortestPath.push_back(activeState);
+	}
+
+	*/
+
+	/*!
+	 *
+	 */
+	static void findShortestPath(storm::storage::SparseMatrix<T> const& transMat, storm::storage::BitVector& subSysStates, storm::storage::BitVector& terminalStates, storm::storage::BitVector& allowedStates, std::vector<std::pair<uint_fast64_t, T>>& itDistances, std::vector<uint_fast64_t>& shortestPath, T& probability) {
+
+		//Do Dijksta to find the shortest path from init states to all states
+		std::vector<std::pair<uint_fast64_t, T>> distances;
+		doDijkstraSearch(transMat, subSysStates, terminalStates, allowedStates, itDistances, distances);
+
+		// Then get the shortest of them
+		extractShortestPath(subSysStates, terminalStates, distances, itDistances, shortestPath, probability,false);
+	}
+
+	/*!
+	 * Only initialized distances vector!
+	 */
+	static void extractShortestPath(storm::storage::BitVector& subSysStates, storm::storage::BitVector& terminalStates, std::vector<std::pair<uint_fast64_t, T>>& distances, std::vector<std::pair<uint_fast64_t, T>>& itDistances, std::vector<uint_fast64_t>& shortestPath, T& probability, bool stopAtFirstTarget = true, bool itSearch = false) {
+		//Find terminal state of best distance
+		uint_fast64_t bestIndex = 0;
+		T bestValue = (T) 0;
+
+		for(storm::storage::BitVector::constIndexIterator term = terminalStates.begin(); term != terminalStates.end(); ++term) {
+
+			//the terminal state might not have been found if it is in a system of forbidden states
+			if(distances[*term].second != -1 && distances[*term].second > bestValue){
+				bestIndex = *term;
+				bestValue = distances[*term].second;
+
+				//if set, stop since the first target that is not null was the only one found
+				if(stopAtFirstTarget) break;
+			}
+		}
+
+		if(!itSearch) {
+			// it's a subSys->subSys search. So target states are terminals and subsys states
+			for(auto term = subSysStates.begin(); term != subSysStates.end(); ++term) {
+
+				//the terminal state might not have been found if it is in a system of forbidden states
+				if(distances[*term].second != -1 && distances[*term].second > bestValue){
+					bestIndex = *term;
+					bestValue = distances[*term].second;
+
+					//if set, stop since the first target that is not null was the only one found
+					if(stopAtFirstTarget) break;
+				}
+			}
+		}
+
+		//safety test: is the best terminal state viable?
+		if(distances[bestIndex].second == (T) -1){
+			shortestPath.push_back(bestIndex);
+			probability = (T) 0;
+			LOG4CPLUS_DEBUG(logger, "Terminal state not viable!");
+			return;
+		}
+
+		// save the probability to reach the state via the shortest path
+		probability = distances[bestIndex].second;
+
+		//Reconstruct shortest path. Notice that the last state of the path might be an initState.
+
+		//There might be a chain of terminal states with 1.0 transitions at the end of the path
+		while(terminalStates.get(distances[bestIndex].first)) {
+			bestIndex = distances[bestIndex].first;
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Found best state: " << bestIndex);
+		LOG4CPLUS_DEBUG(logger, "Value: " << bestValue);
+
+		shortestPath.push_back(bestIndex);
+		bestIndex = distances[bestIndex].first;
+		while(!subSysStates.get(bestIndex)) {
+			shortestPath.push_back(bestIndex);
+			bestIndex = distances[bestIndex].first;
+		}
+		shortestPath.push_back(bestIndex);
+
+		//At last compensate for the distance between init and source state
+		probability = itSearch ? probability : probability / itDistances[bestIndex].second;
+
+		LOG4CPLUS_DEBUG(logger, "Found path: " << shortestPath);
+	}
+
+private:
+
+	template <typename Type>
+	struct transition {
+		uint_fast64_t source;
+		Type prob;
+		uint_fast64_t target;
+	};
+
+	template <class CompareType>
+	class CompareTransitions {
+	public:
+		bool operator()(transition<CompareType>& t1, transition<CompareType>& t2) {
+			return t1.prob < t2.prob;
+		}
+	};
+
+public:
+
+	/*!
+	 *
+	 */
+	static storm::models::Dtmc<T> computeCriticalSubsystem(storm::models::Dtmc<T>& model, storm::property::prctl::AbstractStateFormula<T> const& stateFormula) {
+
+		//-------------------------------------------------------------
+		// 1. Strip and handle formulas
+		//-------------------------------------------------------------
+
+#ifdef BENCHMARK
+		LOG4CPLUS_INFO(logger, "Formula: " << stateFormula.toString());
+#endif
+		LOG4CPLUS_INFO(logger, "Start finding critical subsystem.");
+
+		// make model checker
+		// TODO: Implement and use generic Model Checker factory.
+		storm::modelchecker::prctl::SparseDtmcPrctlModelChecker<T> modelCheck {model, new storm::solver::GmmxxLinearEquationSolver<T>()};
+
+		// init bit vector to contain the subsystem
+		storm::storage::BitVector subSys(model.getNumberOfStates());
+
+		storm::property::prctl::AbstractPathFormula<T> const* pathFormulaPtr;
+		T bound = 0;
+
+		// Strip bound operator
+		storm::property::prctl::ProbabilisticBoundOperator<T> const* boundOperator = dynamic_cast<storm::property::prctl::ProbabilisticBoundOperator<T> const*>(&stateFormula);
+
+		if(boundOperator == nullptr){
+			LOG4CPLUS_ERROR(logger, "No path bound operator at formula root.");
+			return model.getSubDtmc(subSys);
+		}
+		bound = boundOperator->getBound();
+
+		storm::property::prctl::AbstractPathFormula<T> const& abstractPathFormula = boundOperator->getPathFormula();
+		pathFormulaPtr = &abstractPathFormula;
+
+		// get "init" labeled states
+		storm::storage::BitVector initStates = model.getLabeledStates("init");
+
+		//get real prob for formula
+		logger.getAppender("mainFileAppender")->setThreshold(log4cplus::WARN_LOG_LEVEL);
+		std::vector<T> trueProbs = pathFormulaPtr->check(modelCheck, false);
+		logger.getAppender("mainFileAppender")->setThreshold(log4cplus::INFO_LOG_LEVEL);
+
+		T trueProb = 0;
+		for(auto index : initStates) {
+			trueProb += trueProbs[index];
+		}
+		trueProb /= initStates.getNumberOfSetBits();
+		//std::cout << "True Prob: " << trueProb << std::endl;
+
+		// get allowed and target states
+		storm::storage::BitVector allowedStates;
+		storm::storage::BitVector targetStates;
+
+		storm::property::prctl::Eventually<T> const* eventually = dynamic_cast<storm::property::prctl::Eventually<T> const*>(pathFormulaPtr);
+		storm::property::prctl::Globally<T> const* globally = dynamic_cast<storm::property::prctl::Globally<T> const*>(pathFormulaPtr);
+		storm::property::prctl::Until<T> const* until = dynamic_cast<storm::property::prctl::Until<T> const*>(pathFormulaPtr);
+		if(eventually != nullptr) {
+			targetStates = eventually->getChild().check(modelCheck);
+			allowedStates = storm::storage::BitVector(targetStates.getSize(), true);
+		}
+		else if(globally != nullptr){
+			//eventually reaching a state without property visiting only states with property
+			allowedStates = globally->getChild().check(modelCheck);
+			targetStates = storm::storage::BitVector(allowedStates);
+			targetStates.complement();
+		}
+		else if(until != nullptr) {
+			allowedStates = until->getLeft().check(modelCheck);
+			targetStates = until->getRight().check(modelCheck);
+		}
+		else {
+			LOG4CPLUS_ERROR(logger, "Strange path formula. Can't decipher.");
+			return model.getSubDtmc(subSys);
+		}
+
+		//-------------------------------------------------------------
+		// 2. Precomputations for heuristics
+		//-------------------------------------------------------------
+
+		// estimate the path count using the models state count as well as the probability bound
+		uint_fast8_t const minPrec = 10;
+		uint_fast64_t const stateCount = model.getNumberOfStates();
+		uint_fast64_t const stateEstimate = ((T) stateCount) * bound;
+
+		//since this only has a good effect on big models -> use only if model has at least 10^5 states
+		uint_fast64_t precision = stateEstimate > 100000 ? stateEstimate/1000 : minPrec;
+
+
+
+		//-------------------------------------------------------------
+		// 3. Subsystem generation
+		//-------------------------------------------------------------
+
+		// Search from init to target states until the shortest path for each target state is reached.
+		std::vector<uint_fast64_t> shortestPath;
+		std::vector<T> subSysProbs;
+		T pathProb = 0;
+		T subSysProb = 0;
+		uint_fast64_t pathCount = 0;
+		uint_fast64_t mcCount = 0;
+
+		// First test if there are init states that are also target states.
+		// If so the init state represents a subsystem with probability mass 1.
+		// -> return it
+		if((initStates & targetStates).getNumberOfSetBits() != 0) {
+			subSys.set((initStates & targetStates).getSetIndicesList().front(), true);
+
+			LOG4CPLUS_INFO(logger, "Critical subsystem found.");
+			LOG4CPLUS_INFO(logger, "Paths needed: " << pathCount);
+			LOG4CPLUS_INFO(logger, "State count of critical subsystem: " << subSys.getNumberOfSetBits());
+			LOG4CPLUS_INFO(logger, "Prob: " << 1);
+			LOG4CPLUS_INFO(logger, "Model checks: " << mcCount);
+
+			return model.getSubDtmc(subSys);
+		}
+
+
+		// Then compute the shortest paths from init states to all target states
+		std::vector<std::pair<uint_fast64_t, T>> initTargetDistances;
+		computeShortestDistances(model.getTransitionMatrix(), initStates, targetStates, allowedStates, initTargetDistances);
+
+		pathProb = 0;
+		extractShortestPath(initStates, targetStates, initTargetDistances, initTargetDistances, shortestPath, pathProb, false, true);
+
+		// push states of found path into subsystem
+		for(auto index : shortestPath) {
+			subSys.set(index, true);
+		}
+		pathCount++;
+
+		// Get estimate (upper bound) of new sub system probability
+		// That is: prob(target) * cost(path) * (mean(prob(inits))/prob(source))
+		subSysProb += trueProbs[shortestPath[0]] * pathProb * trueProb / trueProbs[shortestPath.back()];
+
+		//find new nodes until the system becomes critical.
+		while(true) {
+			shortestPath.clear();
+			pathProb = 0;
+			findShortestPath(model.getTransitionMatrix(), subSys, targetStates, allowedStates, initTargetDistances, shortestPath, pathProb);
+			//doBackwardsSearch(*model->getTransitionMatrix(), *initStates, subSys, targetStates, allowedStates, trueProbs, shortestPath, pathProb);
+
+			pathCount++;
+
+			// merge found states into subsystem
+			for(uint_fast32_t i = 0; i < shortestPath.size(); i++) {
+				subSys.set(shortestPath[i], true);
+			}
+
+			// Get estimate (upper bound) of new sub system probability
+			// That is: prob(target) * cost(path) * (mean(prob(inits))/prob(source))
+			//subSysProb += (trueProbs[shortestPath.back()] == 0 ? 0 : trueProbs[shortestPath[0]] * pathProb * trueProb / trueProbs[shortestPath.back()]);
+			subSysProb += 1*(trueProbs[shortestPath.back()] == 0 ? 1 : trueProbs[shortestPath[0]] * pathProb == 0 ? 1 : pathProb );
+			//std::cout << "Est. prob: " << subSysProb << std::endl;
+			// Do we want to model check?
+			if((pathCount % precision == 0) && subSysProb >= bound) {
+
+				//get probabilities
+				logger.getAppender("mainFileAppender")->setThreshold(log4cplus::WARN_LOG_LEVEL);
+				subSysProbs = modelCheck.checkUntil(allowedStates & subSys, targetStates & subSys, false);
+				logger.getAppender("mainFileAppender")->setThreshold(log4cplus::INFO_LOG_LEVEL);
+
+				//T diff = subSysProb;
+				//std::cout << "Est. prob: " << diff << std::endl;
+
+				// reset sub system prob to correct value
+				subSysProb = 0;
+				for(auto index : initStates) {
+					subSysProb += subSysProbs[index];
+				}
+				subSysProb /= initStates.getNumberOfSetBits();
+
+				mcCount++;
+				//diff -= subSysProb;
+				//std::cout << "Real prob: " << subSysProb << std::endl;
+				//std::cout << "Diff: " << diff << std::endl;
+				//std::cout << "Path count: " << pathCount << std::endl;
+
+				//Are we critical?
+				if(subSysProb >= bound){
+					break;
+				}
+				else if (stateEstimate > 100000){
+					precision = (stateEstimate/1000) - ((stateEstimate/1000) - minPrec)*(subSysProb/bound);
+				}
+			}
+		}
+
+		LOG4CPLUS_INFO(logger, "Critical subsystem found.");
+		LOG4CPLUS_INFO(logger, "Paths needed: " << pathCount);
+		LOG4CPLUS_INFO(logger, "State count of critical subsystem: " << subSys.getNumberOfSetBits());
+		LOG4CPLUS_INFO(logger, "Prob: " << subSysProb);
+		LOG4CPLUS_INFO(logger, "Model checks: " << mcCount);
+
+		return model.getSubDtmc(subSys);
+	}
+
+    };
+
+    } // namespace counterexamples
+} // namespace storm
+
+#endif /* STORM_COUNTEREXAMPLES_PATHBASEDSUBSYSTEMGENERATOR_H_ */

src/modelchecker/prctl/AbstractModelChecker.h

@@ -21,6 +21,7 @@ namespace prctl {
 #include "src/formula/Prctl.h"
 #include "src/storage/BitVector.h"
 #include "src/models/AbstractModel.h"
+#include "src/settings/Settings.h"
 
 #include "log4cplus/logger.h"
 #include "log4cplus/loggingmacros.h"
@@ -66,7 +67,7 @@ public:
 	/*!
 	 * Constructs an AbstractModelChecker with the given model.
 	 */
-	explicit AbstractModelChecker(storm::models::AbstractModel<Type> const& model) : model(model) {
+	explicit AbstractModelChecker(storm::models::AbstractModel<Type> const& model) : model(model){
 		// Intentionally left empty.
 	}
 	/*!
@@ -157,6 +158,7 @@ public:
 			std::cout << "Error during computation: " << e.what() << "Skipping property." << std::endl;
 			LOG4CPLUS_ERROR(logger, "Error during computation: " << e.what() << "Skipping property.");
 		}
+
 		std::cout << std::endl << "-------------------------------------------" << std::endl;
 	}
 

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@@ -74,6 +74,7 @@ public:
 		return formula.check(*this, false);
 	}
 
+
 	/*!
 	 * Checks the given formula that is a bounded-until formula.
 	 *
@@ -86,14 +87,29 @@ public:
 	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type> checkBoundedUntil(storm::property::prctl::BoundedUntil<Type> const& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the bounded until-formula.
-		storm::storage::BitVector leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector rightStates = formula.getRight().check(*this);
+		return this->checkBoundedUntil(formula.getLeft().check(*this), formula.getRight().check(*this), formula.getBound(), qualitative);
+	}
+
+
+	/*!
+	 * Computes the probability to satisfy phi until psi inside a given bound for each state in the model.
+	 *
+	 * @param phiStates A bit vector indicating which states satisfy phi.
+     * @param psiStates A bit vector indicating which states satisfy psi.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type> checkBoundedUntil(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, bool qualitative) const {
+
         std::vector<Type> result(this->getModel().getNumberOfStates());
         
         // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the
         // further analysis.
-        storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), leftStates, rightStates, true, formula.getBound());
+        storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), phiStates, psiStates, true, stepBound);
         LOG4CPLUS_INFO(logger, "Found " << statesWithProbabilityGreater0.getNumberOfSetBits() << " 'maybe' states.");
         
         // Check if we already know the result (i.e. probability 0) for all initial states and
@@ -111,7 +127,7 @@ public:
             storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0);
             
             // Compute the new set of target states in the reduced system.
-            storm::storage::BitVector rightStatesInReducedSystem = statesWithProbabilityGreater0 % rightStates;
+            storm::storage::BitVector rightStatesInReducedSystem = statesWithProbabilityGreater0 % psiStates;
             
             // Make all rows absorbing that satisfy the second sub-formula.
             submatrix.makeRowsAbsorbing(rightStatesInReducedSystem);
@@ -122,7 +138,7 @@ public:
             
             // Perform the matrix vector multiplication as often as required by the formula bound.
             if (linearEquationSolver != nullptr) {
-                this->linearEquationSolver->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound());
+                this->linearEquationSolver->performMatrixVectorMultiplication(submatrix, subresult, nullptr, stepBound);
             } else {
                 throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available.";
             }
@@ -219,6 +235,7 @@ public:
 		return result;
 	}
 
+
 	/*!
 	 * Check the given formula that is an until formula.
 	 *
@@ -231,13 +248,26 @@ public:
 	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type> checkUntil(storm::property::prctl::Until<Type> const& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
-		storm::storage::BitVector leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector rightStates = formula.getRight().check(*this);
+		return this->checkUntil(formula.getLeft().check(*this), formula.getRight().check(*this), qualitative);
+	}
+
+	/*!
+	 * Computes the  probability to satisfy phi until psi for each state in the model.
+	 *
+	 * @param phiStates A bit vector indicating which states satisfy phi.
+     * @param psiStates A bit vector indicating which states satisfy psi.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type> checkUntil(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
 
-		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
+		// We need to identify the states which have to be taken out of the matrix, i.e.
 		// all states that have probability 0 and 1 of satisfying the until-formula.
-        std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->getModel(), leftStates, rightStates);
+        std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->getModel(), phiStates, psiStates);
 		storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
 		storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
 

src/models/AbstractModel.h

@@ -66,8 +66,8 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * @param transitionMatrix The matrix representing the transitions in the model.
 		 * @param stateLabeling The labeling that assigns a set of atomic
 		 * propositions to each state.
-		 * @param stateRewardVector The reward values associated with the states.
-		 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+		 * @param optionalStateRewardVector The reward values associated with the states.
+		 * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model.
          * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
 		 */
 		AbstractModel(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::models::AtomicPropositionsLabeling const& stateLabeling,
@@ -91,8 +91,9 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * @param transitionMatrix The matrix representing the transitions in the model.
 		 * @param stateLabeling The labeling that assigns a set of atomic
 		 * propositions to each state.
-		 * @param stateRewardVector The reward values associated with the states.
-		 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+		 * @param optionalStateRewardVector The reward values associated with the states.
+		 * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model.
+		 * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
 		 */
 		AbstractModel(storm::storage::SparseMatrix<T>&& transitionMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling,
 				boost::optional<std::vector<T>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,

src/models/AtomicPropositionsLabeling.h

@@ -27,7 +27,7 @@ namespace storm {
 namespace models {
 
 /*!
- * This class manages the labeling of the state space with a fixed number of
+ * This class manages the labeling of the state space with a number of
  * atomic propositions.
  */
 class AtomicPropositionsLabeling {
@@ -130,8 +130,9 @@ public:
 
 	/*!
 	 * Registers the name of a proposition.
-	 * Will throw an error if more atomic propositions are added than were originally declared or if an atomic
-     * proposition is registered twice.
+	 * Will throw an error if an atomic proposition is registered twice.
+	 * If more atomic propositions are added than previously declared, the maximum number of propositions
+	 * is increased and the capacity of the singleLabelings vector is matched with the new maximum.
 	 * 
      * @param ap The name of the atomic proposition to add.
 	 * @return The index of the new proposition.
@@ -142,9 +143,8 @@ public:
 			throw storm::exceptions::OutOfRangeException("Atomic Proposition already exists.");
 		}
 		if (apsCurrent >= apCountMax) {
-			LOG4CPLUS_ERROR(logger, "Added more atomic propositions than previously declared.");
-			throw storm::exceptions::OutOfRangeException("Added more atomic propositions than"
-					"previously declared.");
+			apCountMax++;
+			singleLabelings.reserve(apCountMax);
 		}
 		nameToLabelingMap[ap] = apsCurrent;
         singleLabelings.push_back(storm::storage::BitVector(stateCount));
@@ -226,9 +226,9 @@ public:
 	}
 
 	/*!
-	 * Returns the number of atomic propositions managed by this object (set in the initialization).
-     *
-	 * @return The number of atomic propositions.
+	 * Returns the internal index of a given atomic proposition.
+	 *
+	 * @return The index of the atomic proposition.
 	 */
 	uint_fast64_t getIndexOfProposition(std::string const& ap) const {
 		if (!this->containsAtomicProposition(ap)) {
@@ -284,7 +284,23 @@ public:
 	}
 
 	/*!
-	 * Calculates a hash over all values contained in this Sparse Matrix.
+	 * Adds a state to the labeling.
+	 * Since this operation is quite expensive (resizing of all BitVectors containing the labeling), it should
+	 * only be used in special cases to add one or two states.
+	 * If you want to build a new AtomicPropositionlabeling:
+	 *   - Count the number of states you need.
+	 *   - Then add the labelings using addAtomicProposition() and addAtomicPropositionToState().
+	 * Do NOT use this method for this purpose.
+	 */
+	void addState() {
+		for(uint_fast64_t i = 0; i < apsCurrent; i++) {
+			singleLabelings[i].resize(singleLabelings[i].getSize() + 1);
+		}
+		stateCount++;
+	}
+
+	/*!
+	 * Calculates a hash over all values contained in this Atomic Proposition Labeling.
 	 * @return size_t A Hash Value
 	 */
 	std::size_t getHash() const {

src/models/Dtmc.h

@@ -18,6 +18,7 @@
 #include "src/storage/SparseMatrix.h"
 #include "src/exceptions/InvalidArgumentException.h"
 #include "src/settings/Settings.h"
+#include "src/utility/vector.h"
 
 namespace storm {
 
@@ -38,8 +39,9 @@ public:
 	 * @param probabilityMatrix The matrix representing the transitions in the model.
 	 * @param stateLabeling The labeling that assigns a set of atomic
 	 * propositions to each state.
-	 * @param stateRewardVector The reward values associated with the states.
-	 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+	 * @param optionalStateRewardVector The reward values associated with the states.
+	 * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model.
+	 * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
 	 */
 	Dtmc(storm::storage::SparseMatrix<T> const& probabilityMatrix, storm::models::AtomicPropositionsLabeling const& stateLabeling,
 				boost::optional<std::vector<T>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>> const& optionalTransitionRewardMatrix,
@@ -64,8 +66,9 @@ public:
 	 * @param probabilityMatrix The matrix representing the transitions in the model.
 	 * @param stateLabeling The labeling that assigns a set of atomic
 	 * propositions to each state.
-	 * @param stateRewardVector The reward values associated with the states.
-	 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+	 * @param optionalStateRewardVector The reward values associated with the states.
+	 * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model.
+	 * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
 	 */
 	Dtmc(storm::storage::SparseMatrix<T>&& probabilityMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling,
 				boost::optional<std::vector<T>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,
@@ -121,6 +124,182 @@ public:
 		return AbstractDeterministicModel<T>::getHash();
 	}
 
+	/*!
+	 * Generates a sub-Dtmc of this Dtmc induced by the states specified by the bitvector.
+	 * E.g. a Dtmc that is partial isomorph (on the given states) to this one.
+	 * @param subSysStates A BitVector where all states that should be kept are indicated 
+	 *                     by a set bit of the corresponding index.
+	 *                     Waring: If the vector does not have the correct size, it will be resized.
+	 * @return The sub-Dtmc.
+	 */
+	storm::models::Dtmc<T> getSubDtmc(storm::storage::BitVector& subSysStates) {
+
+
+		// Is there any state in the subsystem?
+		if(subSysStates.getNumberOfSetBits() == 0) {
+			LOG4CPLUS_ERROR(logger, "No states in subsystem!");
+			return storm::models::Dtmc<T>(storm::storage::SparseMatrix<T>(0),
+					  	  	  	  	  	  storm::models::AtomicPropositionsLabeling(this->getStateLabeling(), subSysStates),
+					  	  	  	  	  	  boost::optional<std::vector<T>>(),
+					  	  	  	  	  	  boost::optional<storm::storage::SparseMatrix<T>>(),
+					  	  	  	  	  	  boost::optional<std::vector<std::set<uint_fast64_t>>>());
+		}
+
+		// Does the vector have the right size?
+		if(subSysStates.getSize() != this->getNumberOfStates()) {
+			LOG4CPLUS_INFO(logger, "BitVector has wrong size. Resizing it...");
+			subSysStates.resize(this->getNumberOfStates());
+		}
+
+		// Test if it is a proper subsystem of this Dtmc, i.e. if there is at least one state to be left out.
+		if(subSysStates.getNumberOfSetBits() == subSysStates.getSize()) {
+			LOG4CPLUS_INFO(logger, "All states are kept. This is no proper subsystem.");
+			return storm::models::Dtmc<T>(*this);
+		}
+
+		// 1. Get all necessary information from the old transition matrix
+		storm::storage::SparseMatrix<T> const & origMat = this->getTransitionMatrix();
+		uint_fast64_t const stateCount = origMat.getColumnCount();
+
+		// Iterate over all rows. Count the number of all transitions from the old system to be 
+		// transfered to the new one. Also build a mapping from the state number of the old system 
+		// to the state number of the new system.
+		uint_fast64_t subSysTransitionCount = 0;
+		uint_fast64_t row = 0;
+		uint_fast64_t newRow = 0;
+		std::vector<uint_fast64_t> stateMapping;
+		for(auto iter = origMat.begin(); iter != origMat.end(); ++iter) {
+			if(subSysStates.get(row)){
+				for(auto colIter = iter.begin(); colIter != iter.end(); ++colIter) {
+					if(subSysStates.get(colIter.column())) {
+						subSysTransitionCount++;	
+					} 
+				}
+				stateMapping.push_back(newRow);
+				newRow++;
+			} else {
+				stateMapping.push_back((uint_fast64_t) -1);
+			}
+			row++;
+		}
+
+		// 2. Construct transition matrix
+
+		// Take all states indicated by the vector as well as one additional state s_b as target of
+		// all transitions that target a state that is not kept.
+		uint_fast64_t const newStateCount = subSysStates.getNumberOfSetBits() + 1;
+		storm::storage::SparseMatrix<T> newMat(newStateCount);
+
+		// The number of transitions of the new Dtmc is the number of transitions transfered 
+		// from the old one plus one transition for each state to s_b.
+		newMat.initialize(subSysTransitionCount + newStateCount);
+
+		// Now fill the matrix.
+		newRow = 0;
+		row = 0;
+		T rest = 0;
+		for(auto iter = origMat.begin(); iter != origMat.end(); ++iter) {
+			if(subSysStates.get(row)){
+				// Transfer transitions
+				for(auto colIter = iter.begin(); colIter != iter.end(); ++colIter) {
+					if(subSysStates.get(colIter.column())) {
+						newMat.addNextValue(newRow, stateMapping[colIter.column()], colIter.value());
+					} else {
+						rest += colIter.value();
+					}
+				}
+
+				// Insert the transition taking care of the remaining outgoing probability.
+				newMat.addNextValue(newRow, newStateCount - 1, rest);
+				rest = (T) 0;
+
+				newRow++;
+			}
+			row++;
+		}
+
+		// Insert last transition: self loop on s_b
+		newMat.addNextValue(newStateCount - 1, newStateCount - 1, (T) 1);
+
+		newMat.finalize(false);
+
+		// 3. Take care of the labeling.
+		storm::models::AtomicPropositionsLabeling newLabeling = storm::models::AtomicPropositionsLabeling(this->getStateLabeling(), subSysStates);
+		newLabeling.addState();
+		if(!newLabeling.containsAtomicProposition("s_b")) {
+			newLabeling.addAtomicProposition("s_b");
+		}
+		newLabeling.addAtomicPropositionToState("s_b", newStateCount - 1);
+
+		// 4. Handle the optionals
+
+		boost::optional<std::vector<T>> newStateRewards;
+		if(this->hasStateRewards()) {
+
+			// Get the rewards and move the needed values to the front.
+			std::vector<T> newRewards(this->getStateRewardVector());
+			storm::utility::vector::selectVectorValues(newRewards, subSysStates, newRewards);
+
+			// Throw away all values after the last state and set the reward for s_b to 0.
+			newRewards.resize(newStateCount);
+			newRewards[newStateCount - 1] = (T) 0;
+
+			newStateRewards = newRewards;
+		}
+
+		boost::optional<storm::storage::SparseMatrix<T>> newTransitionRewards;
+		if(this->hasTransitionRewards()) {
+
+			storm::storage::SparseMatrix<T> newTransRewards(newStateCount);
+			newTransRewards.initialize(subSysTransitionCount + newStateCount);
+
+			// Copy the rewards for the kept states
+			newRow = 0;
+			row = 0;
+			for(auto iter = this->getTransitionRewardMatrix().begin(); iter != this->getTransitionRewardMatrix().end(); ++iter) {
+				if(subSysStates.get(row)){
+					// Transfer transition rewards
+					for(auto colIter = iter.begin(); colIter != iter.end(); ++colIter) {
+						if(subSysStates.get(colIter.column())) {
+							newTransRewards.addNextValue(newRow, stateMapping[colIter.column()], colIter.value());
+						}
+					}
+
+					// Insert the reward (e.g. 0) for the transition taking care of the remaining outgoing probability.
+					newTransRewards.addNextValue(newRow, newStateCount - 1, (T) 0);
+
+					newRow++;
+				}
+				row++;
+			}
+
+			newTransitionRewards = newTransRewards;
+		}
+
+		boost::optional<std::vector<std::set<uint_fast64_t>>> newChoiceLabels;
+		if(this->hasChoiceLabels()) {
+
+			// Get the choice label sets and move the needed values to the front.
+			std::vector<std::set<uint_fast64_t>> newChoice(this->getChoiceLabeling());
+			storm::utility::vector::selectVectorValues(newChoice, subSysStates, newChoice);
+
+			// Throw away all values after the last state and set the choice label set for s_b as empty.
+			newChoice.resize(newStateCount);
+			newChoice[newStateCount - 1] = std::set<uint_fast64_t>();
+
+			newChoiceLabels = newChoice;
+		}
+
+		// 5. Make Dtmc from its parts and return it
+		return storm::models::Dtmc<T>(newMat, 
+									  newLabeling,
+									  newStateRewards,
+									  newTransitionRewards,
+									  newChoiceLabels
+									  );
+
+	}
+
 private:
 	/*!
 	 *	@brief Perform some sanity checks.

src/parser/PrctlFileParser.cpp

@@ -15,7 +15,8 @@ namespace parser {
 
 std::list<storm::property::prctl::AbstractPrctlFormula<double>*> PrctlFileParser(std::string filename) {
 	// Open file
-	std::ifstream inputFileStream(filename, std::ios::in);
+	std::ifstream inputFileStream;
+	inputFileStream.open(filename, std::ios::in);
 
 	if (!inputFileStream.is_open()) {
 		std::string message = "Error while opening file ";
@@ -24,16 +25,14 @@ std::list<storm::property::prctl::AbstractPrctlFormula<double>*> PrctlFileParser
 
 	std::list<storm::property::prctl::AbstractPrctlFormula<double>*> result;
 
-	while(!inputFileStream.eof()) {
-		std::string line;
-		//The while loop reads the input file line by line
-		while (std::getline(inputFileStream, line)) {
-			PrctlParser parser(line);
-			if (!parser.parsedComment()) {
-				//lines containing comments will be skipped.
-				LOG4CPLUS_INFO(logger, "Parsed formula \"" + line + "\" into \"" + parser.getFormula()->toString() + "\"");
-				result.push_back(parser.getFormula());
-			}
+	std::string line;
+	//The while loop reads the input file line by line
+	while (std::getline(inputFileStream, line)) {
+		PrctlParser parser(line);
+		if (!parser.parsedComment()) {
+			//lines containing comments will be skipped.
+			LOG4CPLUS_INFO(logger, "Parsed formula \"" + line + "\" into \"" + parser.getFormula()->toString() + "\"");
+			result.push_back(parser.getFormula());
 		}
 	}
 

src/storage/BitVector.h

@@ -166,6 +166,7 @@ public:
         uint_fast64_t nextPosition = 0;
         for (auto position : filter) {
             this->set(nextPosition, other.get(position));
+            nextPosition++;
         }
     }
     

src/storage/SparseMatrix.h

@@ -89,7 +89,7 @@ public:
 	typedef T const* ConstValueIterator;
     
     /*!
-	 * A class representing an iterator over a continguous number of rows of the matrix.
+	 * A class representing an iterator over a continuous number of rows of the matrix.
 	 */
 	class ConstIterator {
 	public:

src/storm.cpp

@@ -1,3 +1,4 @@
+
 /*
  *	STORM - a C++ Rebuild of MRMC
  *	
@@ -28,6 +29,7 @@
 #include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
 #include "src/counterexamples/MILPMinimalLabelSetGenerator.h"
 #include "src/counterexamples/SMTMinimalCommandSetGenerator.h"
+#include "src/counterexamples/PathBasedSubsystemGenerator.h"
 #include "src/parser/AutoParser.h"
 #include "src/parser/PrctlParser.h"
 #include "src/utility/ErrorHandling.h"
@@ -250,12 +252,128 @@ void checkPrctlFormulae(storm::modelchecker::prctl::AbstractModelChecker<double>
 		std::list<storm::property::prctl::AbstractPrctlFormula<double>*> formulaList = storm::parser::PrctlFileParser(chosenPrctlFile);
         
         for (auto formula : formulaList) {
-            modelchecker.check(*formula);
+        	modelchecker.check(*formula);
             delete formula;
         }
 	}
 }
 
+/*!
+ * Handles the counterexample generation control.
+ */
+ void generateCounterExample(storm::parser::AutoParser<double> parser) {
+	LOG4CPLUS_INFO(logger, "Starting counterexample generation.");
+	LOG4CPLUS_INFO(logger, "Testing inputs...");
+
+	//First test output directory.
+	std::string outPath = storm::settings::Settings::getInstance()->getOptionByLongName("counterExample").getArgument(1).getValueAsString();
+	if(outPath.back() != '/' && outPath.back() != '\\') {
+		LOG4CPLUS_ERROR(logger, "The output path is not valid.");
+		return;
+	}
+	std::ofstream testFile(outPath + "test.dot");
+	if(testFile.fail()) {
+		LOG4CPLUS_ERROR(logger, "The output path is not valid.");
+		return;
+	}
+	testFile.close();
+	std::remove((outPath + "test.dot").c_str());
+
+ 	//Differentiate between model types.
+	if(parser.getType() != storm::models::DTMC) {
+		LOG4CPLUS_ERROR(logger, "Counterexample generation for the selected model type is not supported.");
+		return;
+	}
+
+	storm::models::Dtmc<double> model = *parser.getModel<storm::models::Dtmc<double>>();
+	LOG4CPLUS_INFO(logger, "Model is a DTMC.");
+
+	// Get specified PRCTL formulas.
+	std::string const chosenPrctlFile = storm::settings::Settings::getInstance()->getOptionByLongName("counterExample").getArgument(0).getValueAsString();
+	LOG4CPLUS_INFO(logger, "Parsing prctl file: " << chosenPrctlFile << ".");
+	std::list<storm::property::prctl::AbstractPrctlFormula<double>*> formulaList = storm::parser::PrctlFileParser(chosenPrctlFile);
+
+	// Test for each formula if a counterexample can be generated for it.
+	if(formulaList.size() == 0) {
+		LOG4CPLUS_ERROR(logger, "No PRCTL formula file specified.");
+		return;
+	}
+
+	// Get prctl file name without the filetype
+	uint_fast64_t first = 0;
+	if(chosenPrctlFile.find('/') != std::string::npos) {
+		first = chosenPrctlFile.find_last_of('/') + 1;
+	} else if(chosenPrctlFile.find('\\') != std::string::npos) {
+		first = chosenPrctlFile.find_last_of('\\') + 1;
+	}
+
+	uint_fast64_t length;
+	if(chosenPrctlFile.find_last_of('.') != std::string::npos && chosenPrctlFile.find_last_of('.') >= first) {
+		length = chosenPrctlFile.find_last_of('.') - first;
+	} else {
+		length = chosenPrctlFile.length() - first;
+	}
+
+	std::string outFileName = chosenPrctlFile.substr(first, length);
+
+	// Test formulas and do generation
+	uint_fast64_t fIndex = 0;
+	for (auto formula : formulaList) {
+
+		// First check if it is a formula type for which a counterexample can be generated.
+		if (dynamic_cast<storm::property::prctl::AbstractStateFormula<double> const*>(formula) == nullptr) {
+			LOG4CPLUS_ERROR(logger, "Unexpected kind of formula. Expected a state formula.");
+			delete formula;
+			continue;
+		}
+
+		storm::property::prctl::AbstractStateFormula<double> const& stateForm = static_cast<storm::property::prctl::AbstractStateFormula<double> const&>(*formula);
+
+		// Do some output
+		std::cout << "Generating counterexample for formula " << fIndex << ":" << std::endl;
+		LOG4CPLUS_INFO(logger, "Generating counterexample for formula " + std::to_string(fIndex) + ": ");
+		std::cout << "\t" << formula->toString() << "\n" << std::endl;
+		LOG4CPLUS_INFO(logger, formula->toString());
+
+		// Now check if the model does not satisfy the formula.
+		// That is if there is at least one initial state of the model that does not.
+
+		// Also raise the logger threshold for the log file, so that the model check infos aren't logged (useless and there are lots of them)
+		// Lower it again after the model check.
+		logger.getAppender("mainFileAppender")->setThreshold(log4cplus::WARN_LOG_LEVEL);
+		storm::storage::BitVector result = stateForm.check(*createPrctlModelChecker(model));
+		logger.getAppender("mainFileAppender")->setThreshold(log4cplus::INFO_LOG_LEVEL);
+
+		if((result & model.getInitialStates()).getNumberOfSetBits() == model.getInitialStates().getNumberOfSetBits()) {
+			std::cout << "Formula is satisfied. Can not generate counterexample.\n\n" << std::endl;
+			LOG4CPLUS_INFO(logger, "Formula is satisfied. Can not generate counterexample.");
+			delete formula;
+			continue;
+		}
+
+		// Generate counterexample
+		storm::models::Dtmc<double> counterExample = storm::counterexamples::PathBasedSubsystemGenerator<double>::computeCriticalSubsystem(*parser.getModel<storm::models::Dtmc<double>>(), stateForm);
+
+		LOG4CPLUS_INFO(logger, "Found counterexample.");
+
+		// Output counterexample
+		// Do standard output
+		std::cout << "Found counterexample with following properties: " << std::endl;
+		counterExample.printModelInformationToStream(std::cout);
+		std::cout << "For full Dtmc see " << outFileName << "_" << fIndex << ".dot at given output path.\n\n" << std::endl;
+
+		// Write the .dot file
+		std::ofstream outFile(outPath + outFileName + "_" + std::to_string(fIndex) + ".dot");
+		if(outFile.good()) {
+			counterExample.writeDotToStream(outFile, true, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, true);
+			outFile.close();
+		}
+
+		fIndex++;
+		delete formula;
+	}
+ }
+
 /*!
  * Main entry point.
  */
@@ -296,40 +414,47 @@ int main(const int argc, const char* argv[]) {
 
 			storm::parser::AutoParser<double> parser(chosenTransitionSystemFile, chosenLabelingFile, chosenStateRewardsFile, chosenTransitionRewardsFile);
 
-            // Determine which engine is to be used to choose the right model checker.
-			LOG4CPLUS_DEBUG(logger, s->getOptionByLongName("matrixLibrary").getArgument(0).getValueAsString());
-
-			// Depending on the model type, the appropriate model checking procedure is chosen.
-            storm::modelchecker::prctl::AbstractModelChecker<double>* modelchecker = nullptr;
-            parser.getModel<storm::models::AbstractModel<double>>()->printModelInformationToStream(std::cout);
-			switch (parser.getType()) {
-			case storm::models::DTMC:
-				LOG4CPLUS_INFO(logger, "Model is a DTMC.");
-                modelchecker = createPrctlModelChecker(*parser.getModel<storm::models::Dtmc<double>>());
-				checkPrctlFormulae(*modelchecker);
-				break;
-			case storm::models::MDP:
-				LOG4CPLUS_INFO(logger, "Model is an MDP.");
-                modelchecker = createPrctlModelChecker(*parser.getModel<storm::models::Mdp<double>>());
-				checkPrctlFormulae(*modelchecker);
-				break;
-			case storm::models::CTMC:
-                LOG4CPLUS_INFO(logger, "Model is a CTMC.");
-                LOG4CPLUS_ERROR(logger, "The selected model type is not supported.");
-                break;
-			case storm::models::CTMDP:
-                LOG4CPLUS_INFO(logger, "Model is a CTMDP.");
-                LOG4CPLUS_ERROR(logger, "The selected model type is not supported.");
-				break;
-			case storm::models::Unknown:
-			default:
-				LOG4CPLUS_ERROR(logger, "The model type could not be determined correctly.");
-				break;
+			//Should there be a counterexample generated in case the formula is not satisfied?
+			if(s->isSet("counterExample")) {
+
+				generateCounterExample(parser);
+			
+			} else {
+				// Determine which engine is to be used to choose the right model checker.
+				LOG4CPLUS_DEBUG(logger, s->getOptionByLongName("matrixLibrary").getArgument(0).getValueAsString());
+
+				// Depending on the model type, the appropriate model checking procedure is chosen.
+				storm::modelchecker::prctl::AbstractModelChecker<double>* modelchecker = nullptr;
+				parser.getModel<storm::models::AbstractModel<double>>()->printModelInformationToStream(std::cout);
+				switch (parser.getType()) {
+				case storm::models::DTMC:
+					LOG4CPLUS_INFO(logger, "Model is a DTMC.");
+					modelchecker = createPrctlModelChecker(*parser.getModel<storm::models::Dtmc<double>>());
+					checkPrctlFormulae(*modelchecker);
+					break;
+				case storm::models::MDP:
+					LOG4CPLUS_INFO(logger, "Model is an MDP.");
+					modelchecker = createPrctlModelChecker(*parser.getModel<storm::models::Mdp<double>>());
+					checkPrctlFormulae(*modelchecker);
+					break;
+				case storm::models::CTMC:
+					LOG4CPLUS_INFO(logger, "Model is a CTMC.");
+					LOG4CPLUS_ERROR(logger, "The selected model type is not supported.");
+					break;
+				case storm::models::CTMDP:
+					LOG4CPLUS_INFO(logger, "Model is a CTMDP.");
+					LOG4CPLUS_ERROR(logger, "The selected model type is not supported.");
+					break;
+				case storm::models::Unknown:
+				default:
+					LOG4CPLUS_ERROR(logger, "The model type could not be determined correctly.");
+					break;
+				}
+
+				if (modelchecker != nullptr) {
+					delete modelchecker;
+				}
 			}
-            
-            if (modelchecker != nullptr) {
-                delete modelchecker;
-            }
 		} else if (s->isSet("symbolic")) {
             // First, we build the model using the given symbolic model description and constant definitions.
             std::string const& programFile = s->getOptionByLongName("symbolic").getArgument(0).getValueAsString();

src/utility/StormOptions.cpp

@@ -1,38 +1,40 @@
 #include "src/utility/StormOptions.h"
 
 bool storm::utility::StormOptions::optionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* settings) -> bool {
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "help", "h", "Shows all available options, arguments and descriptions.").build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "verbose", "v", "Be verbose.").build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "debug", "", "Be very verbose (intended for debugging).").build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "trace", "", "Be extremly verbose (intended for debugging, heavy performance impacts).").build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "logfile", "l", "If specified, the log output will also be written to this file.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("logFileName", "The path and name of the file to write to.").build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "configfile", "c", "If specified, this file will be read and parsed for additional configuration settings.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("configFileName", "The path and name of the file from which to read.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "explicit", "", "Explicit parsing from transition- and labeling files.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("transitionFileName", "The path and name of the file from which to read the transitions.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).addArgument(storm::settings::ArgumentBuilder::createStringArgument("labelingFileName", "The path and name of the file from which to read the labeling.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "symbolic", "", "Parse the given symbolic model file.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("symbolicFileName", "The path and name of the file from which to read the symbolic model.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "prctl", "", "Performs model checking for the PRCTL formulas given in the file.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("prctlFileName", "The file from which to read the PRCTL formulas.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "csl", "", "Performs model checking for the CSL formulas given in the file.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("cslFileName", "The file from which to read the CSL formulas.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "ltl", "", "Performs model checking for the LTL formulas given in the file.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("ltlFileName", "The file from which to read the LTL formulas.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "transitionRewards", "", "If specified, the transition rewards are read from this file and added to the explicit model. Note that this requires an explicit model.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("transitionRewardsFileName", "The file from which to read the rransition rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "stateRewards", "", "If specified, the state rewards are read from this file and added to the explicit model. Note that this requires an explicit model.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("stateRewardsFileName", "The file from which to read the state rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
-    
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "fixDeadlocks", "", "If the model contains deadlock states, setting this option will insert self-loops for these states.").build());
-    
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "help", "h", "Shows all available Options, Arguments and Descriptions").build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "verbose", "v", "Be verbose").build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "debug", "", "Be very verbose (intended for debugging)").build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "trace", "", "Be extremly verbose (intended for debugging, heavy performance impacts)").build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "logfile", "l", "If specified, the log output will also be written to this file").addArgument(storm::settings::ArgumentBuilder::createStringArgument("logFileName", "The path and name of the File to write to").build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "configfile", "c", "If specified, this file will be read and parsed for additional configuration settings").addArgument(storm::settings::ArgumentBuilder::createStringArgument("configFileName", "The path and name of the File to read from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "explicit", "", "Explicit parsing from Transition- and Labeling Files").addArgument(storm::settings::ArgumentBuilder::createStringArgument("transitionFileName", "The path and name of the File to read the transition system from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).addArgument(storm::settings::ArgumentBuilder::createStringArgument("labelingFileName", "The path and name of the File to read the labeling from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "symbolic", "", "Parse the given PRISM File").addArgument(storm::settings::ArgumentBuilder::createStringArgument("prismFileName", "The path and name of the File to read the PRISM Model from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "prctl", "", "Evaluates the PRCTL Formulas given in the File").addArgument(storm::settings::ArgumentBuilder::createStringArgument("prctlFileName", "The path and name of the File to read PRCTL Formulas from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "csl", "", "Evaluates the CSL Formulas given in the File").addArgument(storm::settings::ArgumentBuilder::createStringArgument("cslFileName", "The path and name of the File to read CSL Formulas from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "ltl", "", "Evaluates the LTL Formulas given in the File").addArgument(storm::settings::ArgumentBuilder::createStringArgument("ltlFileName", "The path and name of the File to read LTL Formulas from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "counterExample", "", "Generates a counterexample for the given formulas if not satisfied by the model").addArgument(storm::settings::ArgumentBuilder::createStringArgument("prctlFileName", "The path and name of the File to read the PRCTL Formulas to be used in \n\t\t\t\t\t\t the couterexample generation from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).addArgument(storm::settings::ArgumentBuilder::createStringArgument("outputPath", "The path to the directory to write the generated counterexample files to.").build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "transitionRewards", "", "If specified, the model will have these transition rewards").addArgument(storm::settings::ArgumentBuilder::createStringArgument("transitionRewardsFileName", "The path and name of the File to read the Transition Rewards from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "stateRewards", "", "If specified, the model will have these state rewards").addArgument(storm::settings::ArgumentBuilder::createStringArgument("stateRewardsFileName", "The path and name of the File to read the State Rewards from").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "fixDeadlocks", "", "Insert Self-Loops for States with no outgoing transitions").build());
+
 	std::vector<std::string> matrixLibrarys;
 	matrixLibrarys.push_back("gmm++");
 	matrixLibrarys.push_back("native");
 	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "matrixLibrary", "m", "Sets which matrix library is preferred for numerical operations.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("matrixLibraryName", "The name of a matrix library. Valid values are gmm++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(matrixLibrarys)).setDefaultValueString("gmm++").build()).build());
 
 	return true;
-});
+});