Set up command flow for subsystem generation. Results seem correct on the first look.

Next up: Write the DTMC output for the subsystem generation.

Also:
- fixed one bug in the subsystem generation leading to an infinite loop
- corrected a typo in a comment in SparseMatrix

Former-commit-id: 6014bf2dd2

src/counterexamples/PathBasedSubsystemGenerator.h

@@ -13,6 +13,7 @@
 #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"
@@ -66,13 +67,13 @@ public:
 						distances[*init].second = (T) 1;
 				}
 
-				typename storm::storage::SparseMatrix<T>::ConstRowIterator trans = transMat.begin(*init);
-				for(; trans != transMat.end(*init); ++trans) {
+				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 = trans.row();
+							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));
@@ -90,7 +91,7 @@ public:
 								}
 							}
 
-							distances[trans.column()].first = trans.row();
+							distances[trans.column()].first = *init;
 							distances[trans.column()].second = trans.value();
 
 							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), trans.value()));
@@ -113,7 +114,8 @@ public:
 			// 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
-				for(typename storm::storage::SparseMatrix<T>::ConstRowIterator trans = transMat.begin(activeState.first); trans != transMat.end(activeState.first); ++trans) {
+				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) {
 
@@ -142,7 +144,7 @@ public:
 								}
 							}
 
-							distances[trans.column()].first = trans.row();
+							distances[trans.column()].first = activeState.first;
 							distances[trans.column()].second = distance;
 							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
 						}
@@ -180,19 +182,19 @@ public:
 						continue;
 				}
 
-				typename storm::storage::SparseMatrix<T>::ConstRowIterator trans = transMat.begin(*init);
-				for(; trans != transMat.end(*init); ++trans) {
+				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 = trans.row();
-							distances[trans.column()].second = trans.value() * (itDistances[trans.row()].second == -1 ? 1 : itDistances[trans.row()].second);
+							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[trans.row()].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.
 
@@ -206,8 +208,8 @@ public:
 							}
 
 							//for initialization of subsys -> subsys search use prob (init -> subsys state -> found state) instead of prob(subsys state -> found state)
-							distances[trans.column()].first = trans.row();
-							distances[trans.column()].second = trans.value() * (itDistances[trans.row()].second == -1 ? 1 : itDistances[trans.row()].second);
+							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()));
 						}
@@ -232,7 +234,8 @@ public:
 			// 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
-				for(typename storm::storage::SparseMatrix<T>::ConstRowIterator trans = transMat.begin(activeState.first); trans != transMat.end(activeState.first); ++trans) {
+				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) {
 
@@ -261,7 +264,7 @@ public:
 								}
 							}
 
-							distances[trans.column()].first = trans.row();
+							distances[trans.column()].first = activeState.first;
 							distances[trans.column()].second = distance;
 							activeSet.insert(std::pair<uint_fast64_t, T>(trans.column(), distance));
 						}
@@ -510,7 +513,7 @@ public:
 	/*!
 	 *
 	 */
-	static void computeCriticalSubsystem(storm::models::Dtmc<T> const model, storm::property::prctl::AbstractStateFormula<T> const& stateFormula) {
+	static void computeCriticalSubsystem(storm::models::Dtmc<T> const& model, storm::property::prctl::AbstractStateFormula<T> const& stateFormula) {
 
 		//-------------------------------------------------------------
 		// 1. Strip and handle formulas
@@ -544,7 +547,7 @@ public:
 		pathFormulaPtr = &abstractPathFormula;
 
 		// get "init" labeled states
-		const storm::storage::BitVector initStates = model.getLabeledStates("init");
+		storm::storage::BitVector initStates = model.getLabeledStates("init");
 
 		//get real prob for formula
 		std::vector<T> trueProbs = pathFormulaPtr->check(modelCheck, false);
@@ -608,7 +611,20 @@ public:
 		uint_fast64_t pathCount = 0;
 		uint_fast64_t mcCount = 0;
 
-		// First compute the shortest paths from init states to all target states
+		// 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: " << subSys.toString());
+			std::cout << "Critical subsystem: " << subSys.toString() << std::endl;
+
+			return;
+		}
+
+
+		// 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);
 
@@ -642,7 +658,7 @@ public:
 			// 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 );
+			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) {

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

@@ -257,6 +257,37 @@ void checkPrctlFormulae(storm::modelchecker::prctl::AbstractModelChecker<double>
 	}
 }
 
+/*!
+ * Checks if a counterexample for the given PRCTL formula can be computed.
+ *
+ * @param formula The formula to be examined.
+ * @param model The model for which the formula is to be examined.
+ * @return true iff a counterexample can be computed
+ */
+bool checkForCounterExampleGeneration(storm::property::prctl::AbstractPrctlFormula<double> &formula, storm::models::Dtmc<double> &model) {
+	// First check if it is a formula type for which a counterexample can be generated.
+	if (dynamic_cast<storm::property::prctl::AbstractNoBoundOperator<double> const*>(&formula) != nullptr) {
+		LOG4CPLUS_ERROR(logger, "Can not generate a counterexample without a given bound.");
+		return false;
+	} else if (dynamic_cast<storm::property::prctl::AbstractStateFormula<double> const*>(&formula) != nullptr) {
+		storm::property::prctl::AbstractStateFormula<double> const& stateForm = static_cast<storm::property::prctl::AbstractStateFormula<double> const&>(formula);
+
+		// 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.
+		storm::storage::BitVector result = stateForm.check(*createPrctlModelChecker(model));
+		if((result & model.getInitialStates()).getNumberOfSetBits() < model.getInitialStates().getNumberOfSetBits()) {
+			return true;
+		} else {
+			LOG4CPLUS_ERROR(logger, "Formula is satisfied. Can not generate counterexample.");
+			return false;
+		}
+
+	} else {
+		LOG4CPLUS_ERROR(logger, "Unexpected kind of formula. Expected either a No Bound or a State formula.");
+		return false;
+	}
+}
+
 /*!
  * Main entry point.
  */
@@ -297,40 +328,76 @@ 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("subSys")) {
+
+				//Differentiate between model types.
+				if(parser.getType() == storm::models::DTMC) {
+					LOG4CPLUS_INFO(logger, "Model is a DTMC.");
+
+					if (s->isSet("prctl")) {
+						// Get specified PRCTL formulas.
+						std::string const chosenPrctlFile = s->getOptionByLongName("prctl").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.
+						for (auto formula : formulaList) {
+							if(checkForCounterExampleGeneration(*formula, *parser.getModel<storm::models::Dtmc<double>>())) {
+								//Generate counterexample
+								storm::counterexamples::PathBasedSubsystemGenerator<double>::computeCriticalSubsystem(*parser.getModel<storm::models::Dtmc<double>>(), static_cast<storm::property::prctl::AbstractStateFormula<double> const&>(*formula));
+
+								//Output counterexample
+								//TODO: Write output.
+
+							} else {
+								LOG4CPLUS_INFO(logger, "No counterexample generated for PRCTL formula: " << formula->toString());
+							}
+
+							delete formula;
+						}
+					} else {
+						LOG4CPLUS_ERROR(logger, "No PRCTL formula specified.");
+					}
+				} else {
+					LOG4CPLUS_ERROR(logger, "Counterexample generation for the selected model type is not supported.");
+				}
+			} 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")) {
 			std::string const& programFile = s->getOptionByLongName("symbolic").getArgument(0).getValueAsString();
 			std::string const& constants = s->getOptionByLongName("constants").getArgument(0).getValueAsString();

src/utility/StormOptions.cpp

@@ -9,9 +9,10 @@ bool storm::utility::StormOptions::optionsRegistered = storm::settings::Settings
 	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()).addArgument(storm::settings::ArgumentBuilder::createBooleanArgument("subSys", "generate subsystem as counterexample if formula is not satisfied").setDefaultValueBoolean(false).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", "subSys", "", "Generate subsystem as counterexample if given formula is not satisfied").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());