diff --git a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp index 02c7d783c..67b2dbf1f 100644 --- a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp @@ -26,7 +26,8 @@ namespace storm { namespace helper { template - std::vector computeValidInitialScheduler(uint64_t numberOfMaybeStates, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& b) { + std::vector computeValidInitialSchedulerForUntilProbabilities(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& b) { + uint64_t numberOfMaybeStates = transitionMatrix.getRowGroupCount(); std::vector result(numberOfMaybeStates); storm::storage::BitVector targetStates(numberOfMaybeStates); @@ -106,7 +107,7 @@ namespace storm { std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); // Create the solution vector. - std::vector x(maybeStates.getNonZeroCount(), storm::utility::zero()); + std::vector x(explicitRepresentation.first.getRowGroupCount(), storm::utility::zero()); // Check for requirements of the solver. storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities, dir); @@ -114,7 +115,7 @@ namespace storm { if (!requirements.empty()) { if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { STORM_LOG_DEBUG("Computing valid scheduler hint, because the solver requires it."); - initialScheduler = computeValidInitialScheduler(x.size(), explicitRepresentation.first, explicitRepresentation.second); + initialScheduler = computeValidInitialSchedulerForUntilProbabilities(explicitRepresentation.first, explicitRepresentation.second); requirements.set(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler, false); } @@ -250,6 +251,69 @@ namespace storm { return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), odd, x)); } + template + storm::storage::BitVector computeTargetStatesForReachabilityRewardsFromExplicitRepresentation(storm::storage::SparseMatrix const& transitionMatrix) { + storm::storage::BitVector targetStates(transitionMatrix.getRowGroupCount()); + + // A state is a target state if its row group is empty. + for (uint64_t rowGroup = 0; rowGroup < transitionMatrix.getRowGroupCount(); ++rowGroup) { + if (transitionMatrix.getRowGroupIndices()[rowGroup] == transitionMatrix.getRowGroupIndices()[rowGroup + 1]) { + targetStates.set(rowGroup); + } + } + + return targetStates; + } + + template + std::vector computeValidInitialSchedulerForReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& properMaybeStates, storm::storage::BitVector const& targetStates) { + uint64_t numberOfMaybeStates = transitionMatrix.getRowGroupCount(); + std::vector result(numberOfMaybeStates); + + storm::storage::Scheduler validScheduler(numberOfMaybeStates); + storm::storage::SparseMatrix backwardTransitions = transitionMatrix.transpose(true); + storm::utility::graph::computeSchedulerProb1E(storm::storage::BitVector(numberOfMaybeStates, true), transitionMatrix, backwardTransitions, properMaybeStates, targetStates, validScheduler); + + for (uint64_t state = 0; state < numberOfMaybeStates; ++state) { + if (!targetStates.get(state)) { + result[state] = validScheduler.getChoice(state).getDeterministicChoice(); + } + } + + return result; + } + + template + void eliminateTargetStatesFromExplicitRepresentation(std::pair, std::vector>& explicitRepresentation, std::vector& scheduler, storm::storage::BitVector const& properMaybeStates) { + // Treat the matrix first. + explicitRepresentation.first = explicitRepresentation.first.getSubmatrix(true, properMaybeStates, properMaybeStates); + + // Now eliminate the superfluous entries from the rhs vector and the scheduler. + uint64_t position = 0; + for (auto state : properMaybeStates) { + explicitRepresentation.second[position] = explicitRepresentation.second[state]; + scheduler[position] = scheduler[state]; + position++; + } + + uint64_t numberOfProperMaybeStates = properMaybeStates.getNumberOfSetBits(); + explicitRepresentation.second.resize(numberOfProperMaybeStates); + explicitRepresentation.second.shrink_to_fit(); + scheduler.resize(numberOfProperMaybeStates); + scheduler.shrink_to_fit(); + } + + template + std::vector insertTargetStateValuesIntoExplicitRepresentation(std::vector const& x, storm::storage::BitVector const& properMaybeStates) { + std::vector expandedResult(properMaybeStates.size(), storm::utility::zero()); + + uint64_t position = 0; + for (auto state : properMaybeStates) { + expandedResult[state] = x[position]; + position++; + } + } + template std::unique_ptr HybridMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { @@ -266,7 +330,8 @@ namespace storm { infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, targetStates, storm::utility::graph::performProbGreater0A(model, transitionMatrixBdd, model.getReachableStates(), targetStates)); } infinityStates = !infinityStates && model.getReachableStates(); - storm::dd::Bdd maybeStates = (!targetStates && !infinityStates) && model.getReachableStates(); + storm::dd::Bdd maybeStatesWithTargetStates = !infinityStates && model.getReachableStates(); + storm::dd::Bdd maybeStates = !targetStates && maybeStatesWithTargetStates; STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states."); STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states."); STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); @@ -279,51 +344,85 @@ namespace storm { } else { // If there are maybe states, we need to solve an equation system. if (!maybeStates.isZero()) { + // Check for requirements of the solver this early so we can adapt the maybe states accordingly. + storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards); + bool requireInitialScheduler = false; + if (!requirements.empty()) { + if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) { + STORM_LOG_DEBUG("Computing valid scheduler hint, because the solver requires it."); + requireInitialScheduler = true; + requirements.set(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler, false); + } + STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver."); + } + + // Compute the set of maybe states that we are required to keep in the translation to explicit. + storm::dd::Bdd requiredMaybeStates = requireInitialScheduler ? maybeStatesWithTargetStates : maybeStates; + // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd = maybeStates.createOdd(); + storm::dd::Odd odd = requiredMaybeStates.createOdd(); // 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; + // Start by getting rid of + // (a) transitions from non-maybe states, and + // (b) the choices in the transition matrix that lead to a state that is neither a maybe state + // nor a target state ('infinity choices'). + storm::dd::Add choiceFilterAdd = (transitionMatrixBdd && maybeStatesWithTargetStates.renameVariables(model.getRowVariables(), model.getColumnVariables())).existsAbstract(model.getColumnVariables()).template toAdd(); + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd * choiceFilterAdd; // Then compute the reward vector to use in the computation. storm::dd::Add subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables()); if (!rewardModel.hasStateActionRewards() && !rewardModel.hasTransitionRewards()) { // If the reward model neither has state-action nor transition rewards, we need to multiply // it with the legal nondetermism encodings in each state. - subvector *= transitionMatrixBdd.existsAbstract(model.getColumnVariables()).template toAdd(); + subvector *= choiceFilterAdd; } - // 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); - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. storm::dd::Add stateActionAdd = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd(); std::vector rowGroupSizes = stateActionAdd.sumAbstract(model.getNondeterminismVariables()).toVector(odd); - // Finally cut away all columns targeting non-maybe states. - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - - // Create the solution vector. - std::vector x(maybeStates.getNonZeroCount(), storm::utility::zero()); + // Finally cut away all columns targeting non-maybe states (or non-(maybe or target) states, respectively). + submatrix *= requireInitialScheduler ? maybeStatesWithTargetStates.swapVariables(model.getRowColumnMetaVariablePairs()).template toAdd() : maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); // Translate the symbolic matrix/vector to their explicit representations. std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); - // Check for requirements of the solver. - storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards); - STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver."); - + // Compute a valid initial scheduler if necessary. + boost::optional> initialScheduler; + boost::optional properMaybeStates; + if (requireInitialScheduler) { + // Compute a valid initial scheduler. + storm::storage::BitVector targetStates = computeTargetStatesForReachabilityRewardsFromExplicitRepresentation(explicitRepresentation.first); + properMaybeStates = ~targetStates; + initialScheduler = computeValidInitialSchedulerForReachabilityRewards(explicitRepresentation.first, properMaybeStates.get(), targetStates); + + // Since we needed the transitions to target states to be translated as well for the computation + // of the scheduler, we have to get rid of them now. + eliminateTargetStatesFromExplicitRepresentation(explicitRepresentation, initialScheduler.get(), properMaybeStates.get()); + } + // Now solve the resulting equation system. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); + + // Move the scheduler to the solver. + if (initialScheduler) { + solver->setInitialScheduler(std::move(initialScheduler.get())); + } + + // Create the solution vector. + std::vector x(explicitRepresentation.first.getRowGroupCount(), storm::utility::zero()); + solver->setRequirementsChecked(); solver->solveEquations(dir, x, explicitRepresentation.second); + // If we included the target states in the ODD, we need to expand the result from the solver. + if (requireInitialScheduler) { + x = insertTargetStateValuesIntoExplicitRepresentation(x, properMaybeStates.get()); + } + // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), model.getManager().template getAddZero()), maybeStates, odd, x)); } else { diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index d842705ac..605f108bf 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -519,7 +519,6 @@ namespace storm { } else { storm::storage::BitVector trueStates(transitionMatrix.getRowGroupCount(), true); if (goal.minimize()) { - STORM_LOG_WARN("Results of reward computation may be too low, because of zero-reward loops."); infinityStates = storm::utility::graph::performProb1E(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, trueStates, targetStates); } else { infinityStates = storm::utility::graph::performProb1A(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, trueStates, targetStates); @@ -566,7 +565,7 @@ namespace storm { } // Obtain proper hint information either from the provided hint or from requirements of the solver. - SparseMdpHintType hintInformation = computeHints(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), targetStates, minMaxLinearEquationSolverFactory, selectedChoices); + SparseMdpHintType hintInformation = computeHints(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, ~targetStates, targetStates, minMaxLinearEquationSolverFactory, selectedChoices); // Now compute the results for the maybe states. MaybeStateResult resultForMaybeStates = computeValuesForMaybeStates(goal, submatrix, b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation); diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index ddac7bd2c..a9ef8d0b0 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -623,24 +623,36 @@ namespace storm { } for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { - *targetIt = *sourceIt; - ++sourceIt; - localChoice = 1; - if (choices != nullptr) { - *choiceIt = 0; - } - - for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { - if (filter(*sourceIt, *targetIt)) { - *targetIt = *sourceIt; - if (choices != nullptr) { - *choiceIt = localChoice; + // Only do work if the row group is not empty. + if (*rowGroupingIt != *(rowGroupingIt + 1)) { + *targetIt = *sourceIt; + ++sourceIt; + localChoice = 1; + if (choices != nullptr) { + *choiceIt = 0; + } + + for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { + if (filter(*sourceIt, *targetIt)) { + *targetIt = *sourceIt; + if (choices != nullptr) { + *choiceIt = localChoice; + } } } - } - - if (choices != nullptr) { - ++choiceIt; + + if (choices != nullptr) { + ++choiceIt; + } + } else { + // Compensate for the 'wrong' move forward in the loop header. + --targetIt; + + // Record dummy choice. + if (choices != nullptr) { + *choiceIt = 0; + ++choiceIt; + } } } }); @@ -657,27 +669,39 @@ namespace storm { } for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { - *targetIt = *sourceIt; - ++sourceIt; - localChoice = 1; - if (choices != nullptr) { - *choiceIt = 0; - } - for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { - if (filter(*sourceIt, *targetIt)) { - *targetIt = *sourceIt; - if (choices != nullptr) { - *choiceIt = localChoice; + // Only do work if the row group is not empty. + if (*rowGroupingIt != *(rowGroupingIt + 1)) { + *targetIt = *sourceIt; + ++sourceIt; + localChoice = 1; + if (choices != nullptr) { + *choiceIt = 0; + } + for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { + if (filter(*sourceIt, *targetIt)) { + *targetIt = *sourceIt; + if (choices != nullptr) { + *choiceIt = localChoice; + } } } - } - if (choices != nullptr) { - ++choiceIt; + if (choices != nullptr) { + ++choiceIt; + } + } else { + // Compensate for the 'wrong' move forward in the loop header. + --targetIt; + + // Record dummy choice. + if (choices != nullptr) { + *choiceIt = 0; + ++choiceIt; + } } } #endif } - + /*! * Reduces the given source vector by selecting the smallest element out of each row group. *