Browse Source

Added proper treatment for both upper bound operators to counterexample generators. Added optional statistics output to SAT-based counterexample generator.

Former-commit-id: 5d471c6d00
tempestpy_adaptions
dehnert 11 years ago
parent
commit
e0fbb5cbea
  1. 1
      src/counterexamples/CounterexampleOptions.cpp
  2. 29
      src/counterexamples/MILPMinimalLabelSetGenerator.h
  3. 81
      src/counterexamples/SMTMinimalCommandSetGenerator.h

1
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;
});

29
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<T> const& labeledMdp, VariableInformation const& variableInformation, double probabilityThreshold) {
static uint_fast64_t assertProbabilityGreaterThanThreshold(GRBenv* env, GRBmodel* model, storm::models::Mdp<T> const& labeledMdp, VariableInformation const& variableInformation, double probabilityThreshold, bool strictBound) {
uint_fast64_t numberOfConstraintsCreated = 0;
int error = 0;
int variableIndex = static_cast<int>(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<T> 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<T> 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<uint_fast64_t> getMinimalLabelSet(storm::models::Mdp<T> const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool checkThresholdFeasible = false, bool includeSchedulerCuts = false) {
static std::set<uint_fast64_t> getMinimalLabelSet(storm::models::Mdp<T> 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<uint_fast64_t> usedLabelSet = getMinimalLabelSet(labeledMdp, phiStates, psiStates, bound, true, true);
std::set<uint_fast64_t> 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<std::chrono::milliseconds>(endTime - startTime).count() << "ms." << std::endl;

81
src/counterexamples/SMTMinimalCommandSetGenerator.h

@ -1227,11 +1227,25 @@ namespace storm {
#endif
public:
static std::pair<std::set<uint_fast64_t>, uint_fast64_t > getMinimalCommandSet(storm::ir::Program program, std::string const& constantDefinitionString, storm::models::Mdp<T> const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool checkThresholdFeasible = false) {
static std::set<uint_fast64_t> getMinimalCommandSet(storm::ir::Program program, std::string const& constantDefinitionString, storm::models::Mdp<T> 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<T> subMdp = labeledMdp.restrictChoiceLabels(commandSet);
storm::modelchecker::prctl::SparseMdpPrctlModelChecker<T> modelchecker(subMdp, new storm::solver::GmmxxNondeterministicLinearEquationSolver<T>());
LOG4CPLUS_DEBUG(logger, "Invoking model checker.");
std::vector<T> 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<std::chrono::seconds>(endTime - iterationTimer).count() >= 5) {
std::cout << "Checked " << iterations << " models in " << std::chrono::duration_cast<std::chrono::seconds>(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::seconds>(std::chrono::high_resolution_clock::now() - localClock).count() >= 5) {
std::cout << "Checked " << iterations << " models in " << std::chrono::duration_cast<std::chrono::seconds>(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<std::chrono::milliseconds>(totalSetupTime).count() << "ms" << std::endl;
std::cout << " * time for solving: " << std::chrono::duration_cast<std::chrono::milliseconds>(totalSolverTime).count() << "ms" << std::endl;
std::cout << " * time for checking: " << std::chrono::duration_cast<std::chrono::milliseconds>(totalModelCheckingTime).count() << "ms" << std::endl;
std::cout << " * time for analysis: " << std::chrono::duration_cast<std::chrono::milliseconds>(totalAnalysisTime).count() << "ms" << std::endl;
std::cout << "------------------------------------------" << std::endl;
std::cout << " * total time: " << std::chrono::duration_cast<std::chrono::milliseconds>(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<double>(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<std::chrono::milliseconds>(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<std::chrono::milliseconds>(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;
Loading…
Cancel
Save