From 9498705c952126d1de9a44f263a1cdebbf5c2da1 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 24 Oct 2017 22:13:59 +0200 Subject: [PATCH] more reuse of values in bisimulation-based abstraction refinement --- .../abstraction/GameBasedMdpModelChecker.h | 2 +- .../PartialBisimulationMdpModelChecker.cpp | 196 ++++++++++++------ .../PartialBisimulationMdpModelChecker.h | 13 ++ .../prctl/helper/SymbolicMdpPrctlHelper.cpp | 178 ++++++++-------- .../prctl/helper/SymbolicMdpPrctlHelper.h | 6 +- .../bisimulation/PartialQuotientExtractor.cpp | 4 +- 6 files changed, 246 insertions(+), 153 deletions(-) diff --git a/src/storm/modelchecker/abstraction/GameBasedMdpModelChecker.h b/src/storm/modelchecker/abstraction/GameBasedMdpModelChecker.h index b42fde67f..8ba64de23 100644 --- a/src/storm/modelchecker/abstraction/GameBasedMdpModelChecker.h +++ b/src/storm/modelchecker/abstraction/GameBasedMdpModelChecker.h @@ -9,7 +9,7 @@ #include "storm/storage/SymbolicModelDescription.h" -#include "storm/abstraction/QualitativeResult.h" +#include "storm/abstraction/QualitativeGameResult.h" #include "storm/abstraction/QualitativeGameResultMinMax.h" #include "storm/logic/Bound.h" diff --git a/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.cpp b/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.cpp index 108b57686..bd7f8c100 100644 --- a/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.cpp +++ b/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.cpp @@ -10,6 +10,7 @@ #include "storm/modelchecker/results/SymbolicQuantitativeCheckResult.h" #include "storm/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h" #include "storm/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h" +#include "storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h" #include "storm/logic/FragmentSpecification.h" @@ -93,24 +94,30 @@ namespace storm { // Obtain lower and upper bounds from the partial quotient. std::pair, std::unique_ptr> bounds = computeBoundsPartialQuotient(*quotient, checkTask); - // If either of the two bounds does not exist, the answer can be derived from the existing bounds. - if (bounds.first == nullptr || bounds.second == nullptr) { - STORM_LOG_ASSERT(bounds.first || bounds.second, "Expected at least one bound."); - quotient->printModelInformationToStream(std::cout); - if (bounds.first) { - return std::move(bounds.first); - } else { - return std::move(bounds.second); + bool converged = false; + if (!bounds.first && !bounds.second) { + STORM_LOG_TRACE("Did not compute any bounds, skipping convergence check."); + } else { + // If either of the two bounds does not exist, the answer can be derived from the existing bounds. + if (bounds.first == nullptr || bounds.second == nullptr) { + STORM_LOG_ASSERT(bounds.first || bounds.second, "Expected at least one bound."); + STORM_LOG_TRACE("Obtained result on partial quotient."); + quotient->printModelInformationToStream(std::cout); + if (bounds.first) { + return std::move(bounds.first); + } else { + return std::move(bounds.second); + } } - } - // Check whether the bounds are sufficiently close. - bool converged = checkBoundsSufficientlyClose(bounds); - if (converged) { - result = getAverageOfBounds(bounds); - } else { - printBoundsInformation(bounds); + // Check whether the bounds are sufficiently close. + converged = checkBoundsSufficientlyClose(bounds); + if (converged) { + result = getAverageOfBounds(bounds); + } + } + if (!converged) { STORM_LOG_TRACE("Performing bisimulation step."); bisimulation.compute(10); } @@ -236,6 +243,95 @@ namespace storm { return result; } + template + storm::abstraction::QualitativeMdpResultMinMax::DdType> PartialBisimulationMdpModelChecker::computeQualitativeResult(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask, storm::dd::Bdd const& constraintStates, storm::dd::Bdd const& targetStates) { + + STORM_LOG_DEBUG("Computing qualitative solution."); + storm::abstraction::QualitativeMdpResultMinMax result; + + auto start = std::chrono::high_resolution_clock::now(); + bool isRewardFormula = checkTask.getFormula().isEventuallyFormula() && checkTask.getFormula().asEventuallyFormula().getContext() == storm::logic::FormulaContext::Reward; + storm::dd::Bdd transitionMatrixBdd = quotient.getTransitionMatrix().notZero(); + if (isRewardFormula) { + result.min.second = storm::utility::graph::performProb1E(quotient, transitionMatrixBdd, constraintStates, targetStates, storm::utility::graph::performProbGreater0E(quotient, transitionMatrixBdd, constraintStates, targetStates)); + result.max.second = storm::utility::graph::performProb1A(quotient, transitionMatrixBdd, targetStates, storm::utility::graph::performProbGreater0A(quotient, transitionMatrixBdd, constraintStates, targetStates)); + } else { + result.min = storm::utility::graph::performProb01Min(quotient, transitionMatrixBdd, constraintStates, targetStates); + result.max = storm::utility::graph::performProb01Max(quotient, transitionMatrixBdd, constraintStates, targetStates); + } + auto end = std::chrono::high_resolution_clock::now(); + + auto timeInMilliseconds = std::chrono::duration_cast(end - start).count(); + STORM_LOG_DEBUG("Computed qualitative solution in " << timeInMilliseconds << "ms."); + + return result; + } + + template + std::unique_ptr PartialBisimulationMdpModelChecker::checkForResult(storm::models::symbolic::Mdp const& quotient, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults, CheckTask const& checkTask) { + std::unique_ptr result; + + bool isRewardFormula = checkTask.getFormula().isEventuallyFormula() && checkTask.getFormula().asEventuallyFormula().getContext() == storm::logic::FormulaContext::Reward; + if (isRewardFormula) { + // In the reachability reward case, we can give an answer if all initial states of the system are infinity + // states in the min result. + if ((quotient.getInitialStates() && !qualitativeResults.min.second) == quotient.getInitialStates()) { + result = std::make_unique>(quotient.getReachableStates(), quotient.getInitialStates(), quotient.getInitialStates().ite(quotient.getManager().getConstant(storm::utility::infinity()), quotient.getManager().template getAddZero())); + } + } else { + // In the reachability probability case, we can give the answer if all initial states are prob1 states + // in the min result or if all initial states are prob0 in the max case. + if ((quotient.getInitialStates() && qualitativeResults.min.second) == quotient.getInitialStates()) { + result = std::make_unique>(quotient.getReachableStates(), quotient.getInitialStates(), quotient.getManager().template getAddOne()); + } else if ((quotient.getInitialStates() && qualitativeResults.max.first) == quotient.getInitialStates()) { + result = std::make_unique>(quotient.getReachableStates(), quotient.getInitialStates(), quotient.getManager().template getAddZero()); + } + } + + return result; + } + + template + bool PartialBisimulationMdpModelChecker::skipQuantitativeSolution(storm::models::symbolic::Mdp const& quotient, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults, CheckTask const& checkTask) { + + bool isRewardFormula = checkTask.getFormula().isEventuallyFormula() && checkTask.getFormula().asEventuallyFormula().getContext() == storm::logic::FormulaContext::Reward; + if (isRewardFormula) { + if ((quotient.getInitialStates() && qualitativeResults.min.second) != (quotient.getInitialStates() && qualitativeResults.max.second)) { + return true; + } + } else { + if ((quotient.getInitialStates() && qualitativeResults.min.first) != (quotient.getInitialStates() && qualitativeResults.max.first)) { + return true; + } else if ((quotient.getInitialStates() && qualitativeResults.min.second) != (quotient.getInitialStates() && qualitativeResults.max.second)) { + return true; + } + } + return false; + } + + template + std::pair, std::unique_ptr> PartialBisimulationMdpModelChecker::computeQuantitativeResult(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask, storm::dd::Bdd const& constraintStates, storm::dd::Bdd const& targetStates, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults) { + + std::pair, std::unique_ptr> result; + + bool isRewardFormula = checkTask.getFormula().isEventuallyFormula() && checkTask.getFormula().asEventuallyFormula().getContext() == storm::logic::FormulaContext::Reward; + if (isRewardFormula) { + storm::dd::Bdd maybeMin = qualitativeResults.min.second && quotient.getReachableStates(); + result.first = storm::modelchecker::helper::SymbolicMdpPrctlHelper::computeReachabilityRewards(storm::OptimizationDirection::Minimize, quotient, quotient.getTransitionMatrix(), quotient.getTransitionMatrix().notZero(), checkTask.isRewardModelSet() ? quotient.getRewardModel(checkTask.getRewardModel()) : quotient.getRewardModel(""), maybeMin, targetStates, !qualitativeResults.min.second && quotient.getReachableStates(), storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory(), quotient.getManager().template getAddZero()); + + storm::dd::Bdd maybeMax = qualitativeResults.max.second && quotient.getReachableStates(); + result.second = storm::modelchecker::helper::SymbolicMdpPrctlHelper::computeReachabilityRewards(storm::OptimizationDirection::Maximize, quotient, quotient.getTransitionMatrix(), quotient.getTransitionMatrix().notZero(), checkTask.isRewardModelSet() ? quotient.getRewardModel(checkTask.getRewardModel()) : quotient.getRewardModel(""), maybeMin, targetStates, !qualitativeResults.min.second && quotient.getReachableStates(), storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory(), maybeMax.ite(result.first->asSymbolicQuantitativeCheckResult().getValueVector(), quotient.getManager().template getAddZero())); + } else { + storm::dd::Bdd maybeMin = !(qualitativeResults.min.first || qualitativeResults.min.second) && quotient.getReachableStates(); + result.first = storm::modelchecker::helper::SymbolicMdpPrctlHelper::computeUntilProbabilities(storm::OptimizationDirection::Minimize, quotient, quotient.getTransitionMatrix(), maybeMin, qualitativeResults.min.second, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory(), quotient.getManager().template getAddZero()); + + storm::dd::Bdd maybeMax = !(qualitativeResults.max.first || qualitativeResults.max.second) && quotient.getReachableStates(); + result.second = storm::modelchecker::helper::SymbolicMdpPrctlHelper::computeUntilProbabilities(storm::OptimizationDirection::Maximize, quotient, quotient.getTransitionMatrix(), maybeMax, qualitativeResults.max.second, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory(), maybeMax.ite(result.first->asSymbolicQuantitativeCheckResult().getValueVector(), quotient.getManager().template getAddZero())); + } + + return result; + } + template std::pair, std::unique_ptr> PartialBisimulationMdpModelChecker::computeBoundsPartialQuotient(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask) { @@ -244,63 +340,31 @@ namespace storm { // We go through two phases. In phase (1) we are solving the qualitative part and in phase (2) the quantitative part. // Preparation: determine the constraint states and the target states of the reachability objective. - bool isRewardFormula = checkTask.getFormula().isEventuallyFormula() && checkTask.getFormula().asEventuallyFormula().getContext() == storm::logic::FormulaContext::Reward; std::pair, storm::dd::Bdd> constraintTargetStates = getConstraintAndTargetStates(quotient, checkTask); // Phase (1): solve qualitatively. - storm::abstraction::QualitativeMdpResultMinMax qualitativeResults; - storm::dd::Bdd transitionMatrixBdd = quotient.getTransitionMatrix().notZero(); - if (isRewardFormula) { - qualitativeResults.min.second = storm::utility::graph::performProb1E(quotient, transitionMatrixBdd, constraintTargetStates.first, constraintTargetStates.second, storm::utility::graph::performProbGreater0E(quotient, transitionMatrixBdd, constraintTargetStates.first, constraintTargetStates.second)); - qualitativeResults.max.second = storm::utility::graph::performProb1A(quotient, transitionMatrixBdd, constraintTargetStates.first, storm::utility::graph::performProbGreater0A(quotient, transitionMatrixBdd, constraintTargetStates.first, constraintTargetStates.second)); - } else { - qualitativeResults.min = storm::utility::graph::performProb01Min(quotient, transitionMatrixBdd, constraintTargetStates.first, constraintTargetStates.second); - qualitativeResults.max = storm::utility::graph::performProb01Max(quotient, transitionMatrixBdd, constraintTargetStates.first, constraintTargetStates.second); + storm::abstraction::QualitativeMdpResultMinMax qualitativeResults = computeQualitativeResult(quotient, checkTask, constraintTargetStates.first, constraintTargetStates.second); + + // Check whether the answer can be given after the qualitative solution. + result.first = checkForResult(quotient, qualitativeResults, checkTask); + if (result.first) { + return result; } + // Check whether we should skip the quantitative solution (for example if there are initial states for which + // the value is already known to be different at this point. + bool doSkipQuantitativeSolution = skipQuantitativeSolution(quotient, qualitativeResults, checkTask); + STORM_LOG_TRACE("" << (doSkipQuantitativeSolution ? "Skipping" : "Not skipping") << " quantitative solution."); -// CheckTask newCheckTask(checkTask); -// -// if (!checkTask.isBoundSet() || storm::logic::isLowerBound(checkTask.getBoundComparisonType())) { -// // Check whether we can answer the query (lower bound above bound). -// result.first = computeBoundsPartialQuotient(checker, quotient, rewards, storm::OptimizationDirection::Minimize, newCheckTask); -// if (checkTask.isBoundSet()) { -// bool sufficient = boundsSufficient(quotient, true, result.first->asQuantitativeCheckResult(), checkTask.getBoundComparisonType(), checkTask.getBoundThreshold()); -// if (sufficient) { -// return result; -// } -// } -// -// // Check whether we can answer the query (upper bound below bound). -// result.second = computeBoundsPartialQuotient(checker, quotient, rewards, storm::OptimizationDirection::Maximize, newCheckTask); -// if (checkTask.isBoundSet()) { -// bool sufficient = boundsSufficient(quotient, false, result.second->asQuantitativeCheckResult(), checkTask.getBoundComparisonType(), checkTask.getBoundThreshold()); -// if (sufficient) { -// result.first = nullptr; -// return result; -// } -// } -// } else { -// // Check whether we can answer the query (upper bound below bound). -// result.second = computeBoundsPartialQuotient(checker, quotient, rewards, storm::OptimizationDirection::Maximize, newCheckTask); -// if (checkTask.isBoundSet()) { -// bool sufficient = boundsSufficient(quotient, false, result.second->asQuantitativeCheckResult(), checkTask.getBoundComparisonType(), checkTask.getBoundThreshold()); -// if (sufficient) { -// return result; -// } -// } -// -// // Check whether we can answer the query (lower bound above bound). -// result.first = computeBoundsPartialQuotient(checker, quotient, rewards, storm::OptimizationDirection::Minimize, newCheckTask); -// if (checkTask.isBoundSet()) { -// bool sufficient = boundsSufficient(quotient, true, result.first->asQuantitativeCheckResult(), checkTask.getBoundComparisonType(), checkTask.getBoundThreshold()); -// if (sufficient) { -// result.second = nullptr; -// return result; -// } -// } -// } - + // Phase (2): solve quantitatively. + if (!doSkipQuantitativeSolution) { + result = computeQuantitativeResult(quotient, checkTask, constraintTargetStates.first, constraintTargetStates.second, qualitativeResults); + + storm::modelchecker::SymbolicQualitativeCheckResult initialStateFilter(quotient.getReachableStates(), quotient.getInitialStates()); + result.first->filter(initialStateFilter); + result.second->filter(initialStateFilter); + printBoundsInformation(result); + } return result; } diff --git a/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.h b/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.h index e18c738be..0502870fe 100644 --- a/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.h +++ b/src/storm/modelchecker/abstraction/PartialBisimulationMdpModelChecker.h @@ -28,6 +28,11 @@ namespace storm { } } + namespace abstraction { + template + struct QualitativeMdpResultMinMax; + } + namespace modelchecker { template class SymbolicMdpPrctlModelChecker; @@ -60,6 +65,14 @@ namespace storm { std::unique_ptr getAverageOfBounds(std::pair, std::unique_ptr> const& bounds); void printBoundsInformation(std::pair, std::unique_ptr> const& bounds); + // Methods related to the qualitative solution. + storm::abstraction::QualitativeMdpResultMinMax computeQualitativeResult(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask, storm::dd::Bdd const& constraintStates, storm::dd::Bdd const& targetStates); + std::unique_ptr checkForResult(storm::models::symbolic::Mdp const& quotient, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults, CheckTask const& checkTask); + bool skipQuantitativeSolution(storm::models::symbolic::Mdp const& quotient, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults, CheckTask const& checkTask); + + // Methods related to the quantitative solution. + std::pair, std::unique_ptr> computeQuantitativeResult(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask, storm::dd::Bdd const& constraintStates, storm::dd::Bdd const& targetStates, storm::abstraction::QualitativeMdpResultMinMax const& qualitativeResults); + // Retrieves the constraint and target states of the quotient wrt. to the formula in the check task. std::pair, storm::dd::Bdd> getConstraintAndTargetStates(storm::models::symbolic::Mdp const& quotient, CheckTask const& checkTask); diff --git a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp index f7311559c..ef6cb561b 100644 --- a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp @@ -36,6 +36,52 @@ namespace storm { return result; } + template + std::unique_ptr SymbolicMdpPrctlHelper::computeUntilProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& maybeStates, storm::dd::Bdd const& statesWithProbability1, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues) { + + // Create the matrix and the vector for the equation system. + storm::dd::Add maybeStatesAdd = maybeStates.template toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add prob1StatesAsColumn = statesWithProbability1.template toAdd(); + prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add subvector = submatrix * prob1StatesAsColumn; + subvector = subvector.sumAbstract(model.getColumnVariables()); + + // Finally cut away all columns targeting non-maybe states. + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + + // Now solve the resulting equation system. + std::unique_ptr> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs()); + + // Check requirements of solver. + storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::UntilProbabilities, dir); + boost::optional> initialScheduler; + if (!requirements.empty()) { + if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { + STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it."); + initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability1); + requirements.clearValidInitialScheduler(); + } + requirements.clearBounds(); + STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver."); + } + if (initialScheduler) { + solver->setInitialScheduler(initialScheduler.get()); + } + solver->setBounds(storm::utility::zero(), storm::utility::one()); + solver->setRequirementsChecked(); + + storm::dd::Add result = solver->solveEquations(dir, startValues ? startValues.get() : maybeStatesAdd.getDdManager().template getAddZero(), subvector); + + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability1.template toAdd() + result)); + } + template std::unique_ptr SymbolicMdpPrctlHelper::computeUntilProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues) { // We need to identify the states which have to be taken out of the matrix, i.e. all states that have @@ -48,7 +94,6 @@ namespace storm { } storm::dd::Bdd maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates(); - STORM_LOG_INFO("Preprocessing: " << statesWithProbability01.first.getNonZeroCount() << " states with probability 0, " << statesWithProbability01.second.getNonZeroCount() << " with probability 1 (" << maybeStates.getNonZeroCount() << " states remaining)."); // Check whether we need to compute exact probabilities for some states. @@ -58,47 +103,7 @@ namespace storm { } else { // If there are maybe states, we need to solve an equation system. if (!maybeStates.isZero()) { - // Create the matrix and the vector for the equation system. - storm::dd::Add maybeStatesAdd = maybeStates.template toAdd(); - - // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting - // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - - // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all - // maybe states. - storm::dd::Add prob1StatesAsColumn = statesWithProbability01.second.template toAdd(); - prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); - storm::dd::Add subvector = submatrix * prob1StatesAsColumn; - subvector = subvector.sumAbstract(model.getColumnVariables()); - - // Finally cut away all columns targeting non-maybe states. - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - - // Now solve the resulting equation system. - std::unique_ptr> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs()); - - // Check requirements of solver. - storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::UntilProbabilities, dir); - boost::optional> initialScheduler; - if (!requirements.empty()) { - if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { - STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it."); - initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability01.second); - requirements.clearValidInitialScheduler(); - } - requirements.clearBounds(); - STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver."); - } - if (initialScheduler) { - solver->setInitialScheduler(initialScheduler.get()); - } - solver->setBounds(storm::utility::zero(), storm::utility::one()); - solver->setRequirementsChecked(); - - storm::dd::Add result = solver->solveEquations(dir, startValues ? maybeStates.ite(startValues.get(), model.getManager().template getAddZero()) : model.getManager().template getAddZero(), subvector); - - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.template toAdd() + result)); + return computeUntilProbabilities(dir, model, transitionMatrix, maybeStates, statesWithProbability01.second, linearEquationSolverFactory, startValues ? maybeStates.ite(startValues.get(), model.getManager().template getAddZero()) : model.getManager().template getAddZero()); } else { return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.template toAdd())); } @@ -188,6 +193,52 @@ namespace storm { return std::unique_ptr(new SymbolicQuantitativeCheckResult(model.getReachableStates(), result)); } + template + std::unique_ptr SymbolicMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& transitionMatrixBdd, RewardModelType const& rewardModel, storm::dd::Bdd const& maybeStates, storm::dd::Bdd const& targetStates, storm::dd::Bdd const& infinityStates, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues) { + // Create the matrix and the vector for the equation system. + storm::dd::Add maybeStatesAdd = maybeStates.template toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the state reward vector to use in the computation. + storm::dd::Add subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables()); + + // Since we are cutting away target and infinity states, we need to account for this by giving + // choices the value infinity that have some successor contained in the infinity states. + storm::dd::Bdd choicesWithInfinitySuccessor = (maybeStates && transitionMatrixBdd && infinityStates.swapVariables(model.getRowColumnMetaVariablePairs())).existsAbstract(model.getColumnVariables()); + subvector = choicesWithInfinitySuccessor.ite(model.getManager().template getInfinity(), subvector); + + // Finally cut away all columns targeting non-maybe states. + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + + // Now solve the resulting equation system. + std::unique_ptr> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs()); + + // Check requirements of solver. + storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::ReachabilityRewards, dir); + boost::optional> initialScheduler; + if (!requirements.empty()) { + if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { + STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it."); + initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates); + requirements.clearValidInitialScheduler(); + } + requirements.clearLowerBounds(); + STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver."); + } + if (initialScheduler) { + solver->setInitialScheduler(initialScheduler.get()); + } + solver->setLowerBound(storm::utility::zero()); + solver->setRequirementsChecked(); + + storm::dd::Add result = solver->solveEquations(dir, startValues ? startValues.get() : maybeStatesAdd.getDdManager().template getAddZero(), subvector); + + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), result))); + } + template std::unique_ptr SymbolicMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues) { @@ -210,48 +261,7 @@ namespace storm { // If there are maybe states, we need to solve an equation system. if (!maybeStates.isZero()) { - // Create the matrix and the vector for the equation system. - storm::dd::Add maybeStatesAdd = maybeStates.template toAdd(); - - // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting - // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - - // Then compute the state reward vector to use in the computation. - storm::dd::Add subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables()); - - // Since we are cutting away target and infinity states, we need to account for this by giving - // choices the value infinity that have some successor contained in the infinity states. - storm::dd::Bdd choicesWithInfinitySuccessor = (maybeStates && transitionMatrixBdd && infinityStates.swapVariables(model.getRowColumnMetaVariablePairs())).existsAbstract(model.getColumnVariables()); - subvector = choicesWithInfinitySuccessor.ite(model.getManager().template getInfinity(), subvector); - - // Finally cut away all columns targeting non-maybe states. - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - - // Now solve the resulting equation system. - std::unique_ptr> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs()); - - // Check requirements of solver. - storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::ReachabilityRewards, dir); - boost::optional> initialScheduler; - if (!requirements.empty()) { - if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { - STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it."); - initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates); - requirements.clearValidInitialScheduler(); - } - requirements.clearLowerBounds(); - STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver."); - } - if (initialScheduler) { - solver->setInitialScheduler(initialScheduler.get()); - } - solver->setLowerBound(storm::utility::zero()); - solver->setRequirementsChecked(); - - storm::dd::Add result = solver->solveEquations(dir, startValues ? maybeStates.ite(startValues.get(), model.getManager().template getAddZero()) : model.getManager().template getAddZero(), subvector); - - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), result))); + return computeReachabilityRewards(dir, model, transitionMatrix, transitionMatrixBdd, rewardModel, maybeStates, targetStates, infinityStates, linearEquationSolverFactory, startValues ? maybeStates.ite(startValues.get(), model.getManager().template getAddZero()) : model.getManager().template getAddZero()); } else { return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), model.getManager().template getAddZero()))); } diff --git a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h index 4bbbe11df..21eaf80d6 100644 --- a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h +++ b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h @@ -26,7 +26,9 @@ namespace storm { static std::unique_ptr computeNextProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates); static std::unique_ptr computeUntilProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues = boost::none); - + + static std::unique_ptr computeUntilProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& maybeStates, storm::dd::Bdd const& statesWithProbability1, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues = boost::none); + static std::unique_ptr computeGloballyProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& psiStates, bool qualitative, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); static std::unique_ptr computeCumulativeRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); @@ -34,6 +36,8 @@ namespace storm { static std::unique_ptr computeInstantaneousRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); static std::unique_ptr computeReachabilityRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues = boost::none); + + static std::unique_ptr computeReachabilityRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& transitionMatrixBdd, RewardModelType const& rewardModel, storm::dd::Bdd const& maybeStates, storm::dd::Bdd const& targetStates, storm::dd::Bdd const& infinityStates, storm::solver::SymbolicGeneralMinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional> const& startValues = boost::none); }; } diff --git a/src/storm/storage/dd/bisimulation/PartialQuotientExtractor.cpp b/src/storm/storage/dd/bisimulation/PartialQuotientExtractor.cpp index bfff746b6..9e5e34b54 100644 --- a/src/storm/storage/dd/bisimulation/PartialQuotientExtractor.cpp +++ b/src/storm/storage/dd/bisimulation/PartialQuotientExtractor.cpp @@ -98,7 +98,9 @@ namespace storm { STORM_LOG_TRACE("Quotient transition matrix extracted in " << std::chrono::duration_cast(end - start).count() << "ms."); storm::dd::Bdd quotientTransitionMatrixBdd = quotientTransitionMatrix.notZero(); - storm::dd::Bdd deadlockStates = !quotientTransitionMatrixBdd.existsAbstract(blockPrimeVariableSet) && reachableStates; + std::set nonSourceVariables; + std::set_union(blockPrimeVariableSet.begin(), blockPrimeVariableSet.end(), model.getRowVariables().begin(), model.getRowVariables().end(), std::inserter(nonSourceVariables, nonSourceVariables.begin())); + storm::dd::Bdd deadlockStates = !quotientTransitionMatrixBdd.existsAbstract(nonSourceVariables) && reachableStates; start = std::chrono::high_resolution_clock::now(); std::unordered_map> quotientRewardModels;