From e0fbb5cbea8bf5e500b6627b83370de8ba185e12 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 22 Oct 2013 16:06:53 +0200 Subject: [PATCH] Added proper treatment for both upper bound operators to counterexample generators. Added optional statistics output to SAT-based counterexample generator. Former-commit-id: 5d471c6d0016c1e5c99966489a076fb8809c4baf --- src/counterexamples/CounterexampleOptions.cpp | 1 + .../MILPMinimalLabelSetGenerator.h | 29 ++++--- .../SMTMinimalCommandSetGenerator.h | 81 ++++++++++++++----- 3 files changed, 76 insertions(+), 35 deletions(-) diff --git a/src/counterexamples/CounterexampleOptions.cpp b/src/counterexamples/CounterexampleOptions.cpp index 2fd18d67f..8fb95adfb 100644 --- a/src/counterexamples/CounterexampleOptions.cpp +++ b/src/counterexamples/CounterexampleOptions.cpp @@ -5,5 +5,6 @@ bool CounterexampleOptionsRegistered = storm::settings::Settings::registerNewMod techniques.push_back("sat"); techniques.push_back("milp"); instance->addOption(storm::settings::OptionBuilder("Counterexample", "mincmd", "", "Computes a counterexample for the given symbolic model in terms of a minimal command set.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("propertyFile", "The file containing the properties for which counterexamples are to be generated.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).addArgument(storm::settings::ArgumentBuilder::createStringArgument("method", "Sets which technique is used to derive the counterexample. Must be either \"milp\" or \"sat\".").setDefaultValueString("sat").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(techniques)).build()).build()); + instance->addOption(storm::settings::OptionBuilder("Counterexample", "stats", "s", "Sets whether to display statistics for certain functionalities.").build()); return true; }); diff --git a/src/counterexamples/MILPMinimalLabelSetGenerator.h b/src/counterexamples/MILPMinimalLabelSetGenerator.h index 4c6628bdd..0eed58615 100644 --- a/src/counterexamples/MILPMinimalLabelSetGenerator.h +++ b/src/counterexamples/MILPMinimalLabelSetGenerator.h @@ -182,7 +182,7 @@ namespace storm { int error = error = GRBsetintparam(env, "OutputFlag", 1); // Enable the following line to restrict Gurobi to one thread only. - // error = error = GRBsetintparam(env, "Threads", 1); + error = error = GRBsetintparam(env, "Threads", 1); // Enable the following line to force Gurobi to be as precise about the binary variables as required by the given precision option. error = GRBsetdblparam(env, "IntFeasTol", storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble()); @@ -591,14 +591,15 @@ namespace storm { * @param labeledMdp The labeled MDP. * @param variableInformation A struct with information about the variables of the model. * @param probabilityThreshold The probability that the subsystem must exceed. + * @param strictBound A flag indicating whether the threshold must be exceeded or only matched. * @return The total number of constraints that were created. */ - static uint_fast64_t assertProbabilityGreaterThanThreshold(GRBenv* env, GRBmodel* model, storm::models::Mdp const& labeledMdp, VariableInformation const& variableInformation, double probabilityThreshold) { + static uint_fast64_t assertProbabilityGreaterThanThreshold(GRBenv* env, GRBmodel* model, storm::models::Mdp const& labeledMdp, VariableInformation const& variableInformation, double probabilityThreshold, bool strictBound) { uint_fast64_t numberOfConstraintsCreated = 0; int error = 0; int variableIndex = static_cast(variableInformation.virtualInitialStateVariableIndex); double coefficient = 1.0; - error = GRBaddconstr(model, 1, &variableIndex, &coefficient, GRB_GREATER_EQUAL, probabilityThreshold + 1e-6, nullptr); + error = GRBaddconstr(model, 1, &variableIndex, &coefficient, GRB_GREATER_EQUAL, strictBound ? probabilityThreshold : probabilityThreshold + 1e-6, nullptr); if (error) { LOG4CPLUS_ERROR(logger, "Unable to assert constraint (" << GRBgeterrormsg(env) << ")."); throw storm::exceptions::InvalidStateException() << "Unable to assert constraint (" << GRBgeterrormsg(env) << ")."; @@ -1096,12 +1097,13 @@ namespace storm { * @param choiceInformation The information about the choices in the model. * @param variableInformation A struct with information about the variables of the model. * @param probabilityThreshold The probability threshold the subsystem is required to exceed. + * @param strictBound A flag indicating whether the threshold must be exceeded or only matched. * @param includeSchedulerCuts If set to true, additional constraints are asserted that reduce the set of * possible choices. */ - static void buildConstraintSystem(GRBenv* env, GRBmodel* model, storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& psiStates, StateInformation const& stateInformation, ChoiceInformation const& choiceInformation, VariableInformation const& variableInformation, double probabilityThreshold, bool includeSchedulerCuts = false) { + static void buildConstraintSystem(GRBenv* env, GRBmodel* model, storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& psiStates, StateInformation const& stateInformation, ChoiceInformation const& choiceInformation, VariableInformation const& variableInformation, double probabilityThreshold, bool strictBound, bool includeSchedulerCuts = false) { // Assert that the reachability probability in the subsystem exceeds the given threshold. - uint_fast64_t numberOfConstraints = assertProbabilityGreaterThanThreshold(env, model, labeledMdp, variableInformation, probabilityThreshold); + uint_fast64_t numberOfConstraints = assertProbabilityGreaterThanThreshold(env, model, labeledMdp, variableInformation, probabilityThreshold, strictBound); LOG4CPLUS_DEBUG(logger, "Asserted that reachability probability exceeds threshold."); // Add constraints that assert the policy takes at most one action in each state. @@ -1263,7 +1265,7 @@ namespace storm { public: - static std::set getMinimalLabelSet(storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool checkThresholdFeasible = false, bool includeSchedulerCuts = false) { + static std::set getMinimalLabelSet(storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool strictBound, bool checkThresholdFeasible = false, bool includeSchedulerCuts = false) { #ifdef STORM_HAVE_GUROBI // (0) Check whether the MDP is indeed labeled. if (!labeledMdp.hasChoiceLabels()) { @@ -1277,8 +1279,8 @@ namespace storm { for (auto state : labeledMdp.getInitialStates()) { maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]); } - if (maximalReachabilityProbability <= probabilityThreshold) { - throw storm::exceptions::InvalidArgumentException() << "Given probability threshold " << probabilityThreshold << " can not be achieved in model with maximal reachability probability of " << maximalReachabilityProbability << "."; + if ((strictBound && maximalReachabilityProbability < probabilityThreshold) || (!strictBound && maximalReachabilityProbability <= probabilityThreshold)) { + throw storm::exceptions::InvalidArgumentException() << "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << "."; } // (2) Identify relevant and problematic states. @@ -1298,7 +1300,7 @@ namespace storm { updateModel(environmentModelPair.first, environmentModelPair.second); // (4.3) Construct constraint system. - buildConstraintSystem(environmentModelPair.first, environmentModelPair.second, labeledMdp, psiStates, stateInformation, choiceInformation, variableInformation, probabilityThreshold, includeSchedulerCuts); + buildConstraintSystem(environmentModelPair.first, environmentModelPair.second, labeledMdp, psiStates, stateInformation, choiceInformation, variableInformation, probabilityThreshold, strictBound, includeSchedulerCuts); // Update model. updateModel(environmentModelPair.first, environmentModelPair.second); @@ -1341,10 +1343,11 @@ namespace storm { LOG4CPLUS_ERROR(logger, "Illegal formula " << probBoundFormula->toString() << " for counterexample generation."); throw storm::exceptions::InvalidPropertyException() << "Illegal formula " << probBoundFormula->toString() << " for counterexample generation."; } - if (probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS) { - LOG4CPLUS_ERROR(logger, "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only strict upper bounds are supported for counterexample generation."); - throw storm::exceptions::InvalidPropertyException() << "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only strict upper bounds are supported for counterexample generation."; + if (probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS && probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS_EQUAL) { + LOG4CPLUS_ERROR(logger, "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only upper bounds are supported for counterexample generation."); + throw storm::exceptions::InvalidPropertyException() << "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only upper bounds are supported for counterexample generation."; } + bool strictBound = probBoundFormula->getComparisonOperator() == storm::property::ComparisonType::LESS; // Now derive the probability threshold we need to exceed as well as the phi and psi states. Simultaneously, check whether the formula is of a valid shape. double bound = probBoundFormula->getBound(); @@ -1372,7 +1375,7 @@ namespace storm { // Delegate the actual computation work to the function of equal name. auto startTime = std::chrono::high_resolution_clock::now(); - std::set usedLabelSet = getMinimalLabelSet(labeledMdp, phiStates, psiStates, bound, true, true); + std::set usedLabelSet = getMinimalLabelSet(labeledMdp, phiStates, psiStates, bound, strictBound, true, true); auto endTime = std::chrono::high_resolution_clock::now(); std::cout << std::endl << "Computed minimal label set of size " << usedLabelSet.size() << " in " << std::chrono::duration_cast(endTime - startTime).count() << "ms." << std::endl; diff --git a/src/counterexamples/SMTMinimalCommandSetGenerator.h b/src/counterexamples/SMTMinimalCommandSetGenerator.h index 97ad9115f..2b9bc8995 100644 --- a/src/counterexamples/SMTMinimalCommandSetGenerator.h +++ b/src/counterexamples/SMTMinimalCommandSetGenerator.h @@ -1227,11 +1227,25 @@ namespace storm { #endif public: - static std::pair, uint_fast64_t > getMinimalCommandSet(storm::ir::Program program, std::string const& constantDefinitionString, storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool checkThresholdFeasible = false) { + static std::set getMinimalCommandSet(storm::ir::Program program, std::string const& constantDefinitionString, storm::models::Mdp const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool strictBound, bool checkThresholdFeasible = false) { #ifdef STORM_HAVE_Z3 - auto startTime = std::chrono::high_resolution_clock::now(); - auto endTime = std::chrono::high_resolution_clock::now(); + // Set up all clocks used for time measurement. + auto totalClock = std::chrono::high_resolution_clock::now(); + auto localClock = std::chrono::high_resolution_clock::now(); + decltype(std::chrono::high_resolution_clock::now() - totalClock) totalTime(0); + + auto setupTimeClock = std::chrono::high_resolution_clock::now(); + decltype(std::chrono::high_resolution_clock::now() - setupTimeClock) totalSetupTime(0); + + auto solverClock = std::chrono::high_resolution_clock::now(); + decltype(std::chrono::high_resolution_clock::now() - solverClock) totalSolverTime(0); + auto modelCheckingClock = std::chrono::high_resolution_clock::now(); + decltype(std::chrono::high_resolution_clock::now() - modelCheckingClock) totalModelCheckingTime(0); + + auto analysisClock = std::chrono::high_resolution_clock::now(); + decltype(std::chrono::high_resolution_clock::now() - analysisClock) totalAnalysisTime(0); + storm::utility::ir::defineUndefinedConstants(program, constantDefinitionString); // (0) Check whether the MDP is indeed labeled. @@ -1246,8 +1260,8 @@ namespace storm { for (auto state : labeledMdp.getInitialStates()) { maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]); } - if (maximalReachabilityProbability <= probabilityThreshold) { - throw storm::exceptions::InvalidArgumentException() << "Given probability threshold " << probabilityThreshold << " can not be achieved in model with maximal reachability probability of " << maximalReachabilityProbability << "."; + if ((strictBound && maximalReachabilityProbability < probabilityThreshold) || (!strictBound && maximalReachabilityProbability <= probabilityThreshold)) { + throw storm::exceptions::InvalidArgumentException() << "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << "."; } // (2) Identify all states and commands that are relevant, because only these need to be considered later. @@ -1276,6 +1290,9 @@ namespace storm { assertSymbolicCuts(program, labeledMdp, variableInformation, relevancyInformation, context, solver); LOG4CPLUS_DEBUG(logger, "Asserted symbolic cuts."); + // As we are done with the setup at this point, stop the clock for the setup time. + totalSetupTime = std::chrono::high_resolution_clock::now() - setupTimeClock; + // (8) Find the smallest set of commands that satisfies all constraints. If the probability of // satisfying phi until psi exceeds the given threshold, the set of labels is minimal and can be returned. // Otherwise, the current solution has to be ruled out and the next smallest solution is retrieved from @@ -1291,16 +1308,20 @@ namespace storm { uint_fast64_t zeroProbabilityCount = 0; do { LOG4CPLUS_DEBUG(logger, "Computing minimal command set."); + solverClock = std::chrono::high_resolution_clock::now(); commandSet = findSmallestCommandSet(context, solver, variableInformation, currentBound); + totalSolverTime += std::chrono::high_resolution_clock::now() - solverClock; LOG4CPLUS_DEBUG(logger, "Computed minimal command set of size " << (commandSet.size() + relevancyInformation.knownLabels.size()) << "."); // Restrict the given MDP to the current set of labels and compute the reachability probability. + modelCheckingClock = std::chrono::high_resolution_clock::now(); commandSet.insert(relevancyInformation.knownLabels.begin(), relevancyInformation.knownLabels.end()); storm::models::Mdp subMdp = labeledMdp.restrictChoiceLabels(commandSet); storm::modelchecker::prctl::SparseMdpPrctlModelChecker modelchecker(subMdp, new storm::solver::GmmxxNondeterministicLinearEquationSolver()); LOG4CPLUS_DEBUG(logger, "Invoking model checker."); std::vector result = modelchecker.checkUntil(false, phiStates, psiStates, false, nullptr); LOG4CPLUS_DEBUG(logger, "Computed model checking results."); + totalModelCheckingTime += std::chrono::high_resolution_clock::now() - modelCheckingClock; // Now determine the maximal reachability probability by checking all initial states. maximalReachabilityProbability = 0; @@ -1308,7 +1329,8 @@ namespace storm { maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]); } - if (maximalReachabilityProbability < probabilityThreshold) { + analysisClock = std::chrono::high_resolution_clock::now(); + if ((strictBound && maximalReachabilityProbability < probabilityThreshold) || (!strictBound && maximalReachabilityProbability <= probabilityThreshold)) { if (maximalReachabilityProbability == 0) { ++zeroProbabilityCount; @@ -1323,22 +1345,36 @@ namespace storm { } else { done = true; } + totalAnalysisTime += (std::chrono::high_resolution_clock::now() - analysisClock); ++iterations; - endTime = std::chrono::high_resolution_clock::now(); - if (std::chrono::duration_cast(endTime - iterationTimer).count() >= 5) { - std::cout << "Checked " << iterations << " models in " << std::chrono::duration_cast(endTime - startTime).count() << "s (out of which " << zeroProbabilityCount << " could not reach the target states). Current command set size is " << commandSet.size() << "." << std::endl; - iterationTimer = std::chrono::high_resolution_clock::now(); + if (std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - localClock).count() >= 5) { + std::cout << "Checked " << iterations << " models in " << std::chrono::duration_cast(totalTime).count() << "s (out of which " << zeroProbabilityCount << " could not reach the target states). Current command set size is " << commandSet.size() << "." << std::endl; + localClock = std::chrono::high_resolution_clock::now(); } } while (!done); - - std::cout << "Checked " << iterations << " models in total out of which " << zeroProbabilityCount << " could not reach the target states." << std::endl; - + + // Compute and emit the time measurements. + totalTime = std::chrono::high_resolution_clock::now() - totalClock; + if (storm::settings::Settings::getInstance()->isSet("stats")) { + std::cout << std::endl; + std::cout << "Time breakdown:" << std::endl; + std::cout << " * time for setup: " << std::chrono::duration_cast(totalSetupTime).count() << "ms" << std::endl; + std::cout << " * time for solving: " << std::chrono::duration_cast(totalSolverTime).count() << "ms" << std::endl; + std::cout << " * time for checking: " << std::chrono::duration_cast(totalModelCheckingTime).count() << "ms" << std::endl; + std::cout << " * time for analysis: " << std::chrono::duration_cast(totalAnalysisTime).count() << "ms" << std::endl; + std::cout << "------------------------------------------" << std::endl; + std::cout << " * total time: " << std::chrono::duration_cast(totalTime).count() << "ms" << std::endl; + std::cout << std::endl; + std::cout << "Other:" << std::endl; + std::cout << " * number of models checked: " << iterations << std::endl; + std::cout << " * number of models that could not reach a target state: " << zeroProbabilityCount << " (" << 100 * static_cast(zeroProbabilityCount)/iterations << "%)" << std::endl << std::endl; + } + // (9) Return the resulting command set after undefining the constants. storm::utility::ir::undefineUndefinedConstants(program); - return std::make_pair(commandSet, iterations); - + return commandSet; #else throw storm::exceptions::NotImplementedException() << "This functionality is unavailable since StoRM has been compiled without support for Z3."; #endif @@ -1353,10 +1389,11 @@ namespace storm { LOG4CPLUS_ERROR(logger, "Illegal formula " << probBoundFormula->toString() << " for counterexample generation."); throw storm::exceptions::InvalidPropertyException() << "Illegal formula " << probBoundFormula->toString() << " for counterexample generation."; } - if (probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS) { - LOG4CPLUS_ERROR(logger, "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only strict upper bounds are supported for counterexample generation."); - throw storm::exceptions::InvalidPropertyException() << "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only strict upper bounds are supported for counterexample generation."; + if (probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS && probBoundFormula->getComparisonOperator() != storm::property::ComparisonType::LESS_EQUAL) { + LOG4CPLUS_ERROR(logger, "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only upper bounds are supported for counterexample generation."); + throw storm::exceptions::InvalidPropertyException() << "Illegal comparison operator in formula " << probBoundFormula->toString() << ". Only upper bounds are supported for counterexample generation."; } + bool strictBound = probBoundFormula->getComparisonOperator() == storm::property::ComparisonType::LESS; // Now derive the probability threshold we need to exceed as well as the phi and psi states. Simultaneously, check whether the formula is of a valid shape. double bound = probBoundFormula->getBound(); @@ -1384,13 +1421,13 @@ namespace storm { // Delegate the actual computation work to the function of equal name. auto startTime = std::chrono::high_resolution_clock::now(); - auto labelSetIterationPair = getMinimalCommandSet(program, constantDefinitionString, labeledMdp, phiStates, psiStates, bound, true); + auto labelSet = getMinimalCommandSet(program, constantDefinitionString, labeledMdp, phiStates, psiStates, bound, strictBound, true); auto endTime = std::chrono::high_resolution_clock::now(); - std::cout << std::endl << "Computed minimal label set of size " << labelSetIterationPair.first.size() << " in " << std::chrono::duration_cast(endTime - startTime).count() << "ms (" << labelSetIterationPair.second << " models tested)." << std::endl; + std::cout << std::endl << "Computed minimal label set of size " << labelSet.size() << " in " << std::chrono::duration_cast(endTime - startTime).count() << "ms." << std::endl; - std::cout << "Resulting program:" << std::endl; + std::cout << "Resulting program:" << std::endl << std::endl; storm::ir::Program restrictedProgram(program); - restrictedProgram.restrictCommands(labelSetIterationPair.first); + restrictedProgram.restrictCommands(labelSet); std::cout << restrictedProgram.toString() << std::endl; std::cout << std::endl << "-------------------------------------------" << std::endl;