From 52d729b1c751f2cde643a0ee893944decbe3d5ef Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 20 Sep 2017 14:09:21 +0200 Subject: [PATCH] upper bounds computation for reachability rewards in sparse MDPs --- .../helper/DsMpiUpperRewardBoundsComputer.cpp | 299 ++++++++++++++++++ .../helper/DsMpiUpperRewardBoundsComputer.h | 98 ++++++ .../prctl/helper/SparseDtmcPrctlHelper.cpp | 142 +-------- .../prctl/helper/SparseMdpPrctlHelper.cpp | 118 +++++-- src/storm/solver/AbstractEquationSolver.cpp | 158 +++++++++ src/storm/solver/AbstractEquationSolver.h | 107 ++++++- .../IterativeMinMaxLinearEquationSolver.cpp | 29 +- src/storm/solver/LinearEquationSolver.cpp | 116 ------- src/storm/solver/LinearEquationSolver.h | 88 +----- .../solver/MinMaxLinearEquationSolver.cpp | 38 +-- src/storm/solver/MinMaxLinearEquationSolver.h | 43 +-- .../solver/NativeLinearEquationSolver.cpp | 2 +- .../ConsecutiveUint64DynamicPriorityQueue.h | 1 + 13 files changed, 756 insertions(+), 483 deletions(-) create mode 100644 src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp create mode 100644 src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h create mode 100644 src/storm/solver/AbstractEquationSolver.cpp diff --git a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp new file mode 100644 index 000000000..03b233e70 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp @@ -0,0 +1,299 @@ +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" + +#include "storm-config.h" + +#include "storm/storage/BitVector.h" +#include "storm/storage/ConsecutiveUint64DynamicPriorityQueue.h" +#include "storm/storage/SparseMatrix.h" + +#include "storm/storage/sparse/StateType.h" + +#include "storm/utility/macros.h" +#include "storm/utility/constants.h" +#include "storm/utility/ConstantsComparator.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + DsMpiDtmcUpperRewardBoundsComputer::DsMpiDtmcUpperRewardBoundsComputer(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) : transitionMatrix(transitionMatrix), originalRewards(rewards), originalOneStepTargetProbabilities(oneStepTargetProbabilities), backwardTransitions(transitionMatrix.transpose()), p(transitionMatrix.getRowGroupCount()), w(transitionMatrix.getRowGroupCount()), rewards(rewards), targetProbabilities(oneStepTargetProbabilities) { + // Intentionally left empty. + } + + template + std::vector DsMpiDtmcUpperRewardBoundsComputer::computeUpperBounds() { + sweep(); + ValueType lambda = computeLambda(); + STORM_LOG_TRACE("DS-MPI computed lambda as " << lambda << "."); + + // Finally compute the upper bounds for the states. + std::vector result(transitionMatrix.getRowGroupCount()); + auto one = storm::utility::one(); + for (storm::storage::sparse::state_type state = 0; state < result.size(); ++state) { + result[state] = w[state] + (one - p[state]) * lambda; + } + +#ifndef NDEBUG + ValueType max = storm::utility::zero(); + uint64_t nonZeroCount = 0; + for (auto const& e : result) { + if (!storm::utility::isZero(e)) { + ++nonZeroCount; + max = std::max(max, e); + } + } + STORM_LOG_TRACE("DS-MPI computed " << nonZeroCount << " non-zero upper bounds and a maximal bound of " << max << "."); +#endif + return result; + } + + template + ValueType DsMpiDtmcUpperRewardBoundsComputer::computeLambda() const { + ValueType lambda = storm::utility::zero(); + for (storm::storage::sparse::state_type state = 0; state < transitionMatrix.getRowGroupCount(); ++state) { + lambda = std::max(lambda, computeLambdaForChoice(state)); + } + return lambda; + } + + template + ValueType DsMpiDtmcUpperRewardBoundsComputer::computeLambdaForChoice(uint64_t choice) const { + ValueType localLambda = storm::utility::zero(); + uint64_t state = this->getStateForChoice(choice); + + // Check whether condition (I) or (II) applies. + ValueType sum = storm::utility::zero(); + for (auto const& e : transitionMatrix.getRow(choice)) { + sum += e.getValue() * p[e.getColumn()]; + } + sum += originalOneStepTargetProbabilities[choice]; + + if (p[state] < sum) { + STORM_LOG_TRACE("Condition (I) does apply for state " << state << " as " << p[state] << " < " << sum << "."); + // Condition (I) applies. + localLambda = sum - p[state]; + ValueType nominator = originalRewards[choice]; + for (auto const& e : transitionMatrix.getRow(choice)) { + nominator += e.getValue() * w[e.getColumn()]; + } + nominator -= w[state]; + localLambda = nominator / localLambda; + } else { + STORM_LOG_TRACE("Condition (I) does not apply for state " << state << " as " << p[state] << " < " << sum << "."); + // Here, condition (II) automatically applies and as the resulting local lambda is 0, we + // don't need to consider it. + +#ifndef NDEBUG + // Actually check condition (II). + ValueType sum = originalRewards[choice]; + for (auto const& e : transitionMatrix.getRow(choice)) { + sum += e.getValue() * w[e.getColumn()]; + } + STORM_LOG_WARN_COND(w[state] >= sum || storm::utility::ConstantsComparator().isEqual(w[state], sum), "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << sum << "."); +#endif + } + + return localLambda; + } + + template + uint64_t DsMpiDtmcUpperRewardBoundsComputer::getStateForChoice(uint64_t choice) const { + return choice; + } + + template + class DsMpiDtmcPriorityLess { + public: + DsMpiDtmcPriorityLess(DsMpiDtmcUpperRewardBoundsComputer const& dsmpi) : dsmpi(dsmpi) { + // Intentionally left empty. + } + + bool operator()(storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { + ValueType pa = dsmpi.targetProbabilities[a]; + ValueType pb = dsmpi.targetProbabilities[b]; + if (pa < pb) { + return true; + } else if (pa == pb) { + return dsmpi.rewards[a] > dsmpi.rewards[b]; + } + return false; + } + + private: + DsMpiDtmcUpperRewardBoundsComputer const& dsmpi; + }; + + template + void DsMpiDtmcUpperRewardBoundsComputer::sweep() { + // Create a priority queue that allows for easy retrieval of the currently best state. + storm::storage::ConsecutiveUint64DynamicPriorityQueue> queue(transitionMatrix.getRowCount(), DsMpiDtmcPriorityLess(*this)); + + // Keep track of visited states. + storm::storage::BitVector visited(p.size()); + + while (!queue.empty()) { + // Get first entry in queue. + storm::storage::sparse::state_type currentState = queue.popTop(); + + // Mark state as visited. + visited.set(currentState); + + // Set weight and probability for the state. + w[currentState] = rewards[currentState]; + p[currentState] = targetProbabilities[currentState]; + + for (auto const& e : backwardTransitions.getRow(currentState)) { + if (visited.get(e.getColumn())) { + continue; + } + + // Update reward/probability values. + rewards[e.getColumn()] += e.getValue() * w[currentState]; + targetProbabilities[e.getColumn()] += e.getValue() * p[currentState]; + + // Increase priority of element. + queue.increase(e.getColumn()); + } + } + } + + template + DsMpiMdpUpperRewardBoundsComputer::DsMpiMdpUpperRewardBoundsComputer(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) : DsMpiDtmcUpperRewardBoundsComputer(transitionMatrix, rewards, oneStepTargetProbabilities), policy(transitionMatrix.getRowCount()) { + + // Create a mapping from choices to states. + // Also pick a choice in each state that maximizes the target probability and minimizes the reward. + choiceToState.resize(transitionMatrix.getRowCount()); + for (uint64_t state = 0; state < transitionMatrix.getRowGroupCount(); ++state) { + uint64_t choice = transitionMatrix.getRowGroupIndices()[state]; + + boost::optional minReward; + ValueType maxProb = storm::utility::zero(); + + for (uint64_t row = choice, endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) { + choiceToState[row] = state; + + if (this->targetProbabilities[row] > maxProb) { + maxProb = this->targetProbabilities[row]; + minReward = this->rewards[row]; + choice = row; + } else if (this->targetProbabilities[row] == maxProb && (!minReward || minReward.get() > this->rewards[row])) { + minReward = this->rewards[row]; + choice = row; + } + } + + setChoiceInState(state, choice); + } + } + + template + ValueType DsMpiMdpUpperRewardBoundsComputer::computeLambda() const { + ValueType lambda = storm::utility::zero(); + for (storm::storage::sparse::state_type state = 0; state < this->transitionMatrix.getRowGroupCount(); ++state) { + lambda = std::max(lambda, this->computeLambdaForChoice(this->getChoiceInState(state))); + } + return lambda; + } + + template + uint64_t DsMpiMdpUpperRewardBoundsComputer::getStateForChoice(uint64_t choice) const { + return choiceToState[choice]; + } + + template + class DsMpiMdpPriorityLess { + public: + DsMpiMdpPriorityLess(DsMpiMdpUpperRewardBoundsComputer const& dsmpi) : dsmpi(dsmpi) { + // Intentionally left empty. + } + + bool operator()(storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { + uint64_t choiceA = dsmpi.getChoiceInState(a); + uint64_t choiceB = dsmpi.getChoiceInState(b); + + ValueType pa = dsmpi.targetProbabilities[choiceA]; + ValueType pb = dsmpi.targetProbabilities[choiceB]; + if (pa < pb) { + return true; + } else if (pa == pb) { + return dsmpi.rewards[choiceB] > dsmpi.rewards[choiceB]; + } + return false; + } + + private: + DsMpiMdpUpperRewardBoundsComputer const& dsmpi; + }; + + template + void DsMpiMdpUpperRewardBoundsComputer::sweep() { + // Create a priority queue that allows for easy retrieval of the currently best state. + storm::storage::ConsecutiveUint64DynamicPriorityQueue> queue(this->transitionMatrix.getRowGroupCount(), DsMpiMdpPriorityLess(*this)); + + // Keep track of visited states. + storm::storage::BitVector visited(this->transitionMatrix.getRowGroupCount()); + + while (!queue.empty()) { + // Get first entry in queue. + storm::storage::sparse::state_type currentState = queue.popTop(); + + // Mark state as visited. + visited.set(currentState); + + // Set weight and probability for the state. + uint64_t choiceInCurrentState = this->getChoiceInState(currentState); + this->w[currentState] = this->rewards[choiceInCurrentState]; + this->p[currentState] = this->targetProbabilities[choiceInCurrentState]; + + for (auto const& choiceEntry : this->backwardTransitions.getRow(currentState)) { + uint64_t predecessor = this->getStateForChoice(choiceEntry.getColumn()); + if (visited.get(predecessor)) { + continue; + } + + // Update reward/probability values. + this->rewards[choiceEntry.getColumn()] += choiceEntry.getValue() * this->w[currentState]; + this->targetProbabilities[choiceEntry.getColumn()] += choiceEntry.getValue() * this->p[currentState]; + + // If the choice is not the one that is currently taken in the predecessor state, we might need + // to update it. + uint64_t currentChoiceInPredecessor = this->getChoiceInState(predecessor); + if (currentChoiceInPredecessor != choiceEntry.getColumn()) { + // Check whether the updates choice now becomes a better choice in the predecessor state. + ValueType const& newTargetProbability = this->targetProbabilities[choiceEntry.getColumn()]; + ValueType const& newReward = this->rewards[choiceEntry.getColumn()]; + ValueType const& currentTargetProbability = this->targetProbabilities[this->getChoiceInState(predecessor)]; + ValueType const& currentReward = this->rewards[this->getChoiceInState(predecessor)]; + + if (newTargetProbability > currentTargetProbability || (newTargetProbability == currentTargetProbability && newReward < currentReward)) { + setChoiceInState(predecessor, choiceEntry.getColumn()); + } + } + + // Notify the priority of a potential increase of the priority of the element. + queue.increase(predecessor); + } + } + } + + template + uint64_t DsMpiMdpUpperRewardBoundsComputer::getChoiceInState(uint64_t state) const { + return policy[state]; + } + + template + void DsMpiMdpUpperRewardBoundsComputer::setChoiceInState(uint64_t state, uint64_t choice) { + policy[state] = choice; + } + + template class DsMpiDtmcUpperRewardBoundsComputer; + template class DsMpiMdpUpperRewardBoundsComputer; + +#ifdef STORM_HAVE_CARL + template class DsMpiDtmcUpperRewardBoundsComputer; + template class DsMpiMdpUpperRewardBoundsComputer; +#endif + } + } +} diff --git a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h new file mode 100644 index 000000000..1f972672e --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h @@ -0,0 +1,98 @@ +#pragma once + +#include + +namespace storm { + namespace storage { + template + class SparseMatrix; + } + + namespace modelchecker { + namespace helper { + + template + class DsMpiDtmcPriorityLess; + + template + class DsMpiDtmcUpperRewardBoundsComputer { + public: + /*! + * Creates an object that can compute upper bounds on the expected rewards for the provided DTMC. + * + * @param transitionMatrix The matrix defining the transitions of the system without the transitions + * that lead directly to the goal state. + * @param rewards The rewards of each state. + * @param oneStepTargetProbabilities For each state the probability to go to a goal state in one step. + */ + DsMpiDtmcUpperRewardBoundsComputer(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities); + + /*! + * Computes upper bounds on the expected rewards. + */ + std::vector computeUpperBounds(); + + protected: + /*! + * Performs a Dijkstra sweep. + */ + virtual void sweep(); + + /*! + * Computes the lambda used for the estimation. + */ + virtual ValueType computeLambda() const; + + /*! + * Computes the lambda just for the provided choice. + */ + ValueType computeLambdaForChoice(uint64_t choice) const; + + /*! + * Retrieves the state associated with the given choice. + */ + virtual uint64_t getStateForChoice(uint64_t choice) const; + + // References to input data. + storm::storage::SparseMatrix const& transitionMatrix; + std::vector const& originalRewards; + std::vector const& originalOneStepTargetProbabilities; + + // Derived from input data. + storm::storage::SparseMatrix backwardTransitions; + + // Data that the algorithm uses internally. + std::vector p; + std::vector w; + std::vector rewards; + std::vector targetProbabilities; + + friend class DsMpiDtmcPriorityLess; + }; + + template + class DsMpiMdpPriorityLess; + + template + class DsMpiMdpUpperRewardBoundsComputer : public DsMpiDtmcUpperRewardBoundsComputer { + public: + /*! + * Creates an object that can compute upper bounds on the expected rewards for the provided DTMC. + */ + DsMpiMdpUpperRewardBoundsComputer(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities); + + private: + virtual void sweep() override; + virtual ValueType computeLambda() const override; + virtual uint64_t getStateForChoice(uint64_t choice) const override; + uint64_t getChoiceInState(uint64_t state) const; + void setChoiceInState(uint64_t state, uint64_t choice); + + std::vector choiceToState; + std::vector policy; + + friend class DsMpiMdpPriorityLess; + }; + } + } +} diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index c0374fe7d..6a04dfaca 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -14,6 +14,7 @@ #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/modelchecker/hints/ExplicitModelCheckerHint.h" +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" #include "storm/utility/macros.h" #include "storm/utility/ConstantsComparator.h" @@ -220,150 +221,11 @@ namespace storm { targetStates, qualitative, linearEquationSolverFactory, hint); } - template - class DsMpi { - public: - DsMpi(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) : transitionMatrix(transitionMatrix), originalRewards(rewards), originalOneStepTargetProbabilities(oneStepTargetProbabilities), backwardTransitions(transitionMatrix.transpose()), p(transitionMatrix.getRowCount()), w(transitionMatrix.getRowCount()), rewards(rewards), targetProbabilities(oneStepTargetProbabilities) { - // Intentionally left empty. - } - - std::vector computeUpperBounds() { - sweep(); - ValueType lambda = computeLambda(); - STORM_LOG_TRACE("DS-MPI computed lambda as " << lambda << "."); - - // Finally compute the upper bounds for the states. - std::vector result(transitionMatrix.getRowCount()); - auto one = storm::utility::one(); - for (storm::storage::sparse::state_type state = 0; state < result.size(); ++state) { - result[state] = w[state] + (one - p[state]) * lambda; - } - -#ifndef NDEBUG - ValueType max = storm::utility::zero(); - uint64_t nonZeroCount = 0; - for (auto const& e : result) { - if (!storm::utility::isZero(e)) { - ++nonZeroCount; - max = std::max(max, e); - } - } - STORM_LOG_TRACE("DS-MPI computed " << nonZeroCount << " non-zero upper bounds and a maximal bound of " << max << "."); -#endif - return result; - } - - private: - ValueType computeLambda() { - ValueType lambda = storm::utility::zero(); - for (storm::storage::sparse::state_type state = 0; state < targetProbabilities.size(); ++state) { - // Check whether condition (I) or (II) applies. - ValueType sum = storm::utility::zero(); - for (auto const& e : transitionMatrix.getRow(state)) { - sum += e.getValue() * p[e.getColumn()]; - } - sum += originalOneStepTargetProbabilities[state]; - - if (p[state] < sum) { - // Condition (I) applies. - ValueType localLambda = sum - p[state]; - ValueType nominator = originalRewards[state]; - for (auto const& e : transitionMatrix.getRow(state)) { - nominator += e.getValue() * w[e.getColumn()]; - } - nominator -= w[state]; - localLambda = nominator / localLambda; - lambda = std::max(lambda, localLambda); - } else { - // Here, condition (II) automatically applies and as the resulting local lambda is 0, we - // don't need to consider it. - -#ifndef NDEBUG - // Actually check condition (II). - ValueType sum = originalRewards[state]; - for (auto const& e : transitionMatrix.getRow(state)) { - sum += e.getValue() * w[e.getColumn()]; - } - STORM_LOG_WARN_COND(w[state] >= sum || storm::utility::ConstantsComparator().isEqual(w[state], sum), "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << sum << "."); -#endif - } - } - return lambda; - } - - class PriorityLess { - public: - PriorityLess(DsMpi const& dsmpi) : dsmpi(dsmpi) { - // Intentionally left empty. - } - - bool operator()(storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { - ValueType pa = dsmpi.targetProbabilities[a]; - ValueType pb = dsmpi.targetProbabilities[b]; - if (pa < pb) { - return true; - } else if (pa == pb) { - return dsmpi.rewards[a] > dsmpi.rewards[b]; - } - return false; - } - - private: - DsMpi const& dsmpi; - }; - - void sweep() { - // Create a priority queue that allows for easy retrieval of the currently best state. - storm::storage::ConsecutiveUint64DynamicPriorityQueue queue(transitionMatrix.getRowCount(), PriorityLess(*this)); - - storm::storage::BitVector visited(p.size()); - - while (!queue.empty()) { - // Get first entry in queue. - storm::storage::sparse::state_type currentState = queue.popTop(); - - // Mark state as visited. - visited.set(currentState); - - // Set weight and probability for the state. - w[currentState] = rewards[currentState]; - p[currentState] = targetProbabilities[currentState]; - - for (auto const& e : backwardTransitions.getRow(currentState)) { - if (visited.get(e.getColumn())) { - continue; - } - - // Update reward/probability values. - rewards[e.getColumn()] += e.getValue() * w[currentState]; - targetProbabilities[e.getColumn()] += e.getValue() * p[currentState]; - - // Increase priority of element. - queue.increase(e.getColumn()); - } - } - } - - // References to input data. - storm::storage::SparseMatrix const& transitionMatrix; - std::vector const& originalRewards; - std::vector const& originalOneStepTargetProbabilities; - - // Derived from input data. - storm::storage::SparseMatrix backwardTransitions; - - // Data that the algorithm uses internally. - std::vector p; - std::vector w; - std::vector rewards; - std::vector targetProbabilities; - }; - // This function computes an upper bound on the reachability rewards (see Baier et al, CAV'17). template std::vector computeUpperRewardBounds(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) { std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now(); - DsMpi dsmpi(transitionMatrix, rewards, oneStepTargetProbabilities); + DsMpiDtmcUpperRewardBoundsComputer dsmpi(transitionMatrix, rewards, oneStepTargetProbabilities); std::vector bounds = dsmpi.computeUpperBounds(); std::chrono::high_resolution_clock::time_point end = std::chrono::high_resolution_clock::now(); STORM_LOG_TRACE("Computed upper bounds on rewards in " << std::chrono::duration_cast(end - start).count() << "ms."); diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 56a02d7ea..5b36459d7 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -2,9 +2,9 @@ #include - #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/modelchecker/hints/ExplicitModelCheckerHint.h" +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" #include "storm/models/sparse/StandardRewardModel.h" @@ -110,7 +110,7 @@ namespace storm { template struct SparseMdpHintType { - SparseMdpHintType() : eliminateEndComponents(false) { + SparseMdpHintType() : eliminateEndComponents(false), computeUpperBounds(false) { // Intentionally left empty. } @@ -133,10 +133,22 @@ namespace storm { bool hasUpperResultBound() const { return static_cast(upperResultBound); } + + bool hasUpperResultBounds() const { + return static_cast(upperResultBounds); + } ValueType const& getUpperResultBound() const { return upperResultBound.get(); } + + std::vector& getUpperResultBounds() { + return upperResultBounds.get(); + } + + std::vector const& getUpperResultBounds() const { + return upperResultBounds.get(); + } std::vector& getSchedulerHint() { return schedulerHint.get(); @@ -150,11 +162,17 @@ namespace storm { return eliminateEndComponents; } + bool getComputeUpperBounds() { + return computeUpperBounds; + } + boost::optional> schedulerHint; boost::optional> valueHint; boost::optional lowerResultBound; boost::optional upperResultBound; + boost::optional> upperResultBounds; bool eliminateEndComponents; + bool computeUpperBounds; }; template @@ -239,6 +257,10 @@ namespace storm { } else if (type == storm::solver::EquationSystemType::ReachabilityRewards) { requirements.clearLowerBounds(); } + if (requirements.requiresUpperBounds()) { + result.computeUpperBounds = true; + requirements.clearUpperBounds(); + } STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "There are unchecked requirements of the solver."); } else { STORM_LOG_DEBUG("Solver has no requirements."); @@ -260,6 +282,11 @@ namespace storm { if (!result.hasUpperResultBound() && type == storm::solver::EquationSystemType::UntilProbabilities) { result.upperResultBound = storm::utility::one(); } + + // If we received an upper bound, we can drop the requirement to compute one. + if (result.hasUpperResultBound()) { + result.computeUpperBounds = false; + } return result; } @@ -298,6 +325,9 @@ namespace storm { if (hint.hasUpperResultBound()) { solver->setUpperBound(hint.getUpperResultBound()); } + if (hint.hasUpperResultBounds()) { + solver->setUpperBounds(std::move(hint.getUpperResultBounds())); + } if (hint.hasSchedulerHint()) { solver->setInitialScheduler(std::move(hint.getSchedulerHint())); } @@ -309,6 +339,17 @@ namespace storm { // Solve the corresponding system of equations. solver->solveEquations(x, b); +#ifndef NDEBUG + // As a sanity check, make sure our local upper bounds were in fact correct. + if (solver->hasUpperBound(storm::solver::AbstractEquationSolver::BoundType::Local)) { + auto resultIt = x.begin(); + for (auto const& entry : solver->getUpperBounds()) { + STORM_LOG_ASSERT(*resultIt <= entry, "Expecting result value for state " << std::distance(x.begin(), resultIt) << " to be <= " << entry << ", but got " << *resultIt << "."); + ++resultIt; + } + } +#endif + // Create result. MaybeStateResult result(std::move(x)); @@ -485,23 +526,29 @@ namespace storm { }; template - SparseMdpEndComponentInformation eliminateEndComponents(storm::storage::MaximalEndComponentDecomposition const& endComponentDecomposition, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const* sumColumns, storm::storage::BitVector const* selectedChoices, std::vector const* summand, storm::storage::SparseMatrix& submatrix, std::vector& b) { + SparseMdpEndComponentInformation eliminateEndComponents(storm::storage::MaximalEndComponentDecomposition const& endComponentDecomposition, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const* sumColumns, storm::storage::BitVector const* selectedChoices, std::vector const* summand, storm::storage::SparseMatrix& submatrix, std::vector* columnSumVector, std::vector* summandResultVector) { SparseMdpEndComponentInformation result(endComponentDecomposition, maybeStates); // (1) Compute the number of maybe states not in ECs before any other maybe state. std::vector maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices(); - - // Create temporary vector storing possible transitions to ECs. - std::vector> ecValuePairs; - - // (2) Create the parts of the submatrix and vector b that belong to states not contained in ECs. uint64_t numberOfStates = result.numberOfMaybeStatesNotInEc + result.numberOfEc; STORM_LOG_TRACE("Found " << numberOfStates << " states, " << result.numberOfMaybeStatesNotInEc << " not in ECs, " << result.numberOfMaybeStatesInEc << " in ECs and " << result.numberOfEc << " ECs."); + // Prepare builder and vector storage. storm::storage::SparseMatrixBuilder builder(0, numberOfStates, 0, true, true, numberOfStates); - b.resize(numberOfStates); + STORM_LOG_ASSERT((sumColumns && columnSumVector) || (!sumColumns && !columnSumVector), "Expecting a bit vector for which columns to sum iff there is a column sum result vector."); + if (columnSumVector) { + columnSumVector->resize(numberOfStates); + } + STORM_LOG_ASSERT((summand && summandResultVector) || (!summand && !summandResultVector), "Expecting summand iff there is a summand result vector."); + if (summandResultVector) { + summandResultVector->resize(numberOfStates); + } + std::vector> ecValuePairs; + + // (2) Create the parts of the submatrix and vector b that belong to states not contained in ECs. uint64_t currentRow = 0; for (auto state : maybeStates) { if (!result.isStateInEc(state)) { @@ -515,11 +562,11 @@ namespace storm { ecValuePairs.clear(); if (summand) { - b[currentRow] += (*summand)[row]; + (*summandResultVector)[currentRow] += (*summand)[row]; } for (auto const& e : transitionMatrix.getRow(row)) { if (sumColumns && sumColumns->get(e.getColumn())) { - b[currentRow] += e.getValue(); + (*columnSumVector)[currentRow] += e.getValue(); } else if (maybeStates.get(e.getColumn())) { // If the target state of the transition is not contained in an EC, we can just add the entry. if (result.isStateInEc(e.getColumn())) { @@ -563,11 +610,11 @@ namespace storm { ecValuePairs.clear(); if (summand) { - b[currentRow] += (*summand)[row]; + (*summandResultVector)[currentRow] += (*summand)[row]; } for (auto const& e : transitionMatrix.getRow(row)) { if (sumColumns && sumColumns->get(e.getColumn())) { - b[currentRow] += e.getValue(); + (*columnSumVector)[currentRow] += e.getValue(); } else if (maybeStates.get(e.getColumn())) { // If the target state of the transition is not contained in an EC, we can just add the entry. if (result.isStateInEc(e.getColumn())) { @@ -605,7 +652,7 @@ namespace storm { // Only do more work if there are actually end-components. if (!endComponentDecomposition.empty()) { STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " ECs."); - return eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, b); + return eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, &b, nullptr); } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); computeFixedPointSystemUntilProbabilities(transitionMatrix, qualitativeStateSets, submatrix, b); @@ -893,21 +940,27 @@ namespace storm { } template - void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b) { + void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, std::vector* oneStepTargetProbabilities = nullptr) { // Remove rows and columns from the original transition probability matrix for states whose reward values are already known. // If there are infinity states, we additionally have to remove choices of maybeState that lead to infinity. if (qualitativeStateSets.infinityStates.empty()) { submatrix = transitionMatrix.getSubmatrix(true, qualitativeStateSets.maybeStates, qualitativeStateSets.maybeStates, false); b = totalStateRewardVectorGetter(submatrix.getRowCount(), transitionMatrix, qualitativeStateSets.maybeStates); + if (oneStepTargetProbabilities) { + (*oneStepTargetProbabilities) = transitionMatrix.getConstrainedRowGroupSumVector(qualitativeStateSets.maybeStates, targetStates); + } } else { submatrix = transitionMatrix.getSubmatrix(false, *selectedChoices, qualitativeStateSets.maybeStates, false); b = totalStateRewardVectorGetter(transitionMatrix.getRowCount(), transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true)); storm::utility::vector::filterVectorInPlace(b, *selectedChoices); + if (oneStepTargetProbabilities) { + (*oneStepTargetProbabilities) = transitionMatrix.getConstrainedRowSumVector(*selectedChoices, targetStates); + } } } template - boost::optional> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b) { + boost::optional> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, boost::optional>& oneStepTargetProbabilities) { // Start by computing the choices with reward 0, as we only want ECs within this fragment. storm::storage::BitVector zeroRewardChoices(transitionMatrix.getRowCount()); @@ -951,14 +1004,26 @@ namespace storm { // Only do more work if there are actually end-components. if (doDecomposition && !endComponentDecomposition.empty()) { STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " ECs."); - return eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, b); + return eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, oneStepTargetProbabilities ? &targetStates : nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr, &b); } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); - computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, selectedChoices, totalStateRewardVectorGetter, submatrix, b); + computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr); return boost::none; } } + template + void computeUpperRewardBounds(SparseMdpHintType& hintInformation, storm::OptimizationDirection const& direction, storm::storage::SparseMatrix const& submatrix, std::vector const& choiceRewards, std::vector const& oneStepTargetProbabilities) { + + // For the min-case, we use DS-MPI, for the max-case variant 2 of the Baier et al. paper (CAV'17). + if (direction == storm::OptimizationDirection::Minimize) { + DsMpiMdpUpperRewardBoundsComputer dsmpi(submatrix, choiceRewards, oneStepTargetProbabilities); + hintInformation.upperResultBounds = dsmpi.computeUpperBounds(); + } else { + STORM_LOG_ASSERT(false, "Not yet implemented."); + } + } + template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeReachabilityRewardsHelper(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint) { @@ -997,17 +1062,30 @@ namespace storm { // Declare the components of the equation system we will solve. storm::storage::SparseMatrix submatrix; std::vector b; + + // If we need to compute upper bounds on the reward values, we need the one step probabilities + // to a target state. + boost::optional> oneStepTargetProbabilities; + if (hintInformation.getComputeUpperBounds()) { + oneStepTargetProbabilities = std::vector(); + } // If the hint information tells us that we have to eliminate MECs, we do so now. boost::optional> ecInformation; if (hintInformation.getEliminateEndComponents()) { - ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, selectedChoices, totalStateRewardVectorGetter, submatrix, b); + ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities); // Make sure we are not supposed to produce a scheduler if we actually eliminate end components. STORM_LOG_THROW(!ecInformation || !ecInformation.get().eliminatedEndComponents || !produceScheduler, storm::exceptions::NotSupportedException, "Producing schedulers is not supported if end-components need to be eliminated for the solver."); } else { // Otherwise, we compute the standard equations. - computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, selectedChoices, totalStateRewardVectorGetter, submatrix, b); + computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b); + } + + // If we need to compute upper bounds, do so now. + if (hintInformation.getComputeUpperBounds()) { + STORM_LOG_ASSERT(oneStepTargetProbabilities, "Expecting one step target probability vector to be available."); + computeUpperRewardBounds(hintInformation, goal.direction(), submatrix, b, oneStepTargetProbabilities.get()); } // Now compute the results for the maybe states. diff --git a/src/storm/solver/AbstractEquationSolver.cpp b/src/storm/solver/AbstractEquationSolver.cpp new file mode 100644 index 000000000..95e73e82e --- /dev/null +++ b/src/storm/solver/AbstractEquationSolver.cpp @@ -0,0 +1,158 @@ +#include "storm/solver/AbstractEquationSolver.h" + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/adapters/RationalFunctionAdapter.h" + +#include "storm/utility/macros.h" +#include "storm/exceptions/UnmetRequirementException.h" + +namespace storm { + namespace solver { + + template + void AbstractEquationSolver::setTerminationCondition(std::unique_ptr> terminationCondition) { + this->terminationCondition = std::move(terminationCondition); + } + + template + void AbstractEquationSolver::resetTerminationCondition() { + this->terminationCondition = nullptr; + } + + template + bool AbstractEquationSolver::hasCustomTerminationCondition() const { + return static_cast(this->terminationCondition); + } + + template + TerminationCondition const& AbstractEquationSolver::getTerminationCondition() const { + return *terminationCondition; + } + + template + bool AbstractEquationSolver::hasLowerBound(BoundType const& type) const { + if (type == BoundType::Any) { + return static_cast(lowerBound) || static_cast(lowerBounds); + } else if (type == BoundType::Global) { + return static_cast(lowerBound); + } else if (type == BoundType::Local) { + return static_cast(lowerBounds); + } + return false; + } + + template + bool AbstractEquationSolver::hasUpperBound(BoundType const& type) const { + if (type == BoundType::Any) { + return static_cast(upperBound) || static_cast(upperBounds); + } else if (type == BoundType::Global) { + return static_cast(upperBound); + } else if (type == BoundType::Local) { + return static_cast(upperBounds); + } + return false; + } + + template + void AbstractEquationSolver::setLowerBound(ValueType const& value) { + lowerBound = value; + } + + template + void AbstractEquationSolver::setUpperBound(ValueType const& value) { + upperBound = value; + } + + template + void AbstractEquationSolver::setBounds(ValueType const& lower, ValueType const& upper) { + setLowerBound(lower); + setUpperBound(upper); + } + + template + ValueType const& AbstractEquationSolver::getLowerBound() const { + return lowerBound.get(); + } + + template + ValueType const& AbstractEquationSolver::getUpperBound() const { + return upperBound.get(); + } + + template + std::vector const& AbstractEquationSolver::getLowerBounds() const { + return lowerBounds.get(); + } + + template + std::vector const& AbstractEquationSolver::getUpperBounds() const { + return upperBounds.get(); + } + + template + void AbstractEquationSolver::setLowerBounds(std::vector const& values) { + lowerBounds = values; + } + + template + void AbstractEquationSolver::setUpperBounds(std::vector const& values) { + upperBounds = values; + } + + template + void AbstractEquationSolver::setUpperBounds(std::vector&& values) { + upperBounds = std::move(values); + } + + template + void AbstractEquationSolver::setBounds(std::vector const& lower, std::vector const& upper) { + setLowerBounds(lower); + setUpperBounds(upper); + } + + template + void AbstractEquationSolver::createLowerBoundsVector(std::vector& lowerBoundsVector) const { + if (this->hasLowerBound(BoundType::Local)) { + lowerBoundsVector = this->getLowerBounds(); + } else { + STORM_LOG_THROW(this->hasLowerBound(BoundType::Global), storm::exceptions::UnmetRequirementException, "Cannot create lower bounds vector without lower bound."); + for (auto& e : lowerBoundsVector) { + e = this->getLowerBound(); + } + } + } + + template + void AbstractEquationSolver::createUpperBoundsVector(std::unique_ptr>& upperBoundsVector, uint64_t length) const { + STORM_LOG_ASSERT(this->hasUpperBound(), "Expecting upper bound(s)."); + if (!upperBoundsVector) { + if (this->hasUpperBound(BoundType::Local)) { + STORM_LOG_ASSERT(length == this->getUpperBounds().size(), "Mismatching sizes."); + upperBoundsVector = std::make_unique>(this->getUpperBounds()); + } else { + upperBoundsVector = std::make_unique>(length, this->getUpperBound()); + } + } else { + if (this->hasUpperBound(BoundType::Global)) { + for (auto& e : *upperBoundsVector) { + e = this->getUpperBound(); + } + } else { + auto upperBoundsIt = this->getUpperBounds().begin(); + for (auto& e : *upperBoundsVector) { + e = *upperBoundsIt; + ++upperBoundsIt; + } + } + } + } + + template class AbstractEquationSolver; + +#ifdef STORM_HAVE_CARL + template class AbstractEquationSolver; + template class AbstractEquationSolver; +#endif + + } +} diff --git a/src/storm/solver/AbstractEquationSolver.h b/src/storm/solver/AbstractEquationSolver.h index 907e529fa..09ced5229 100644 --- a/src/storm/solver/AbstractEquationSolver.h +++ b/src/storm/solver/AbstractEquationSolver.h @@ -1,9 +1,12 @@ #ifndef STORM_SOLVER_ABSTRACTEQUATIONSOLVER_H_ #define STORM_SOLVER_ABSTRACTEQUATIONSOLVER_H_ -#include "storm/solver/TerminationCondition.h" #include +#include + +#include "storm/solver/TerminationCondition.h" + namespace storm { namespace solver { @@ -16,36 +19,114 @@ namespace storm { * * @param terminationCondition An object that can be queried whether to terminate early or not. */ - void setTerminationCondition(std::unique_ptr> terminationCondition) { - this->terminationCondition = std::move(terminationCondition); - } + void setTerminationCondition(std::unique_ptr> terminationCondition); /*! * Removes a previously set custom termination condition. */ - void resetTerminationCondition() { - this->terminationCondition = nullptr; - } + void resetTerminationCondition(); /*! * Retrieves whether a custom termination condition has been set. */ - bool hasCustomTerminationCondition() const { - return static_cast(this->terminationCondition); - } + bool hasCustomTerminationCondition() const; /*! * Retrieves the custom termination condition (if any was set). * * @return The custom termination condition. */ - TerminationCondition const& getTerminationCondition() const { - return *terminationCondition; - } + TerminationCondition const& getTerminationCondition() const; + + enum class BoundType { + Global, + Local, + Any + }; + + /*! + * Retrieves whether this solver has a lower bound. + */ + bool hasLowerBound(BoundType const& type = BoundType::Any) const; + + /*! + * Retrieves whether this solver has an upper bound. + */ + bool hasUpperBound(BoundType const& type = BoundType::Any) const; + + /*! + * Sets a lower bound for the solution that can potentially be used by the solver. + */ + void setLowerBound(ValueType const& value); + + /*! + * Sets an upper bound for the solution that can potentially be used by the solver. + */ + void setUpperBound(ValueType const& value); + + /*! + * Sets bounds for the solution that can potentially be used by the solver. + */ + void setBounds(ValueType const& lower, ValueType const& upper); + + /*! + * Retrieves the lower bound (if there is any). + */ + ValueType const& getLowerBound() const; + + /*! + * Retrieves the upper bound (if there is any). + */ + ValueType const& getUpperBound() const; + + /*! + * Retrieves a vector containing the lower bounds (if there are any). + */ + std::vector const& getLowerBounds() const; + + /*! + * Retrieves a vector containing the upper bounds (if there are any). + */ + std::vector const& getUpperBounds() const; + + /*! + * Sets lower bounds for the solution that can potentially be used by the solver. + */ + void setLowerBounds(std::vector const& values); + + /*! + * Sets upper bounds for the solution that can potentially be used by the solver. + */ + void setUpperBounds(std::vector const& values); + + /*! + * Sets upper bounds for the solution that can potentially be used by the solver. + */ + void setUpperBounds(std::vector&& values); + + /*! + * Sets bounds for the solution that can potentially be used by the solver. + */ + void setBounds(std::vector const& lower, std::vector const& upper); protected: + void createUpperBoundsVector(std::unique_ptr>& upperBoundsVector, uint64_t length) const; + void createLowerBoundsVector(std::vector& lowerBoundsVector) const; + // A termination condition to be used (can be unset). std::unique_ptr> terminationCondition; + + // A lower bound if one was set. + boost::optional lowerBound; + + // An upper bound if one was set. + boost::optional upperBound; + + // Lower bounds if they were set. + boost::optional> lowerBounds; + + // Lower bounds if they were set. + boost::optional> upperBounds; }; } diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index cc7ab660c..0d37526d9 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -380,37 +380,12 @@ namespace storm { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - if (this->hasInitialScheduler()) { - STORM_LOG_TRACE("Solving initial scheduler hint."); - // Resolve the nondeterminism according to the initial scheduler. - bool convertToEquationSystem = this->linearEquationSolverFactory->getEquationProblemFormat() == LinearEquationSolverProblemFormat::EquationSystem; - storm::storage::SparseMatrix submatrix = this->A->selectRowsFromRowGroups(this->getInitialScheduler(), convertToEquationSystem); - if (convertToEquationSystem) { - submatrix.convertToEquationSystem(); - } - storm::utility::vector::selectVectorValues(*auxiliaryRowGroupVector, this->getInitialScheduler(), this->A->getRowGroupIndices(), b); - - // Solve the resulting equation system. - auto submatrixSolver = this->linearEquationSolverFactory->create(std::move(submatrix)); - submatrixSolver->setCachingEnabled(true); - if (this->lowerBound) { - submatrixSolver->setLowerBound(this->lowerBound.get()); - } - if (this->upperBound) { - submatrixSolver->setUpperBound(this->upperBound.get()); - } - submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector); - } - // Allow aliased multiplications. bool useGaussSeidelMultiplication = this->linEqSolverA->supportsGaussSeidelMultiplication() && settings.getValueIterationMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel; - // Initialize upper bound vector. - for (auto& e : *this->auxiliaryRowGroupVector) { - e = this->getUpperBound(); - } - std::vector* lowerX = &x; + this->createLowerBoundsVector(*lowerX); + this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount()); std::vector* upperX = this->auxiliaryRowGroupVector.get(); std::vector* tmp = nullptr; if (!useGaussSeidelMultiplication) { diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index 22c75794b..8d00d0178 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -132,122 +132,6 @@ namespace storm { cachedRowVector.reset(); } - template - bool LinearEquationSolver::hasLowerBound(BoundType const& type) const { - if (type == BoundType::Any) { - return static_cast(lowerBound) || static_cast(lowerBounds); - } else if (type == BoundType::Global) { - return static_cast(lowerBound); - } else if (type == BoundType::Local) { - return static_cast(lowerBounds); - } - return false; - } - - template - bool LinearEquationSolver::hasUpperBound(BoundType const& type) const { - if (type == BoundType::Any) { - return static_cast(upperBound) || static_cast(upperBounds); - } else if (type == BoundType::Global) { - return static_cast(upperBound); - } else if (type == BoundType::Local) { - return static_cast(upperBounds); - } - return false; - } - - template - void LinearEquationSolver::setLowerBound(ValueType const& value) { - lowerBound = value; - } - - template - void LinearEquationSolver::setUpperBound(ValueType const& value) { - upperBound = value; - } - - template - void LinearEquationSolver::setBounds(ValueType const& lower, ValueType const& upper) { - setLowerBound(lower); - setUpperBound(upper); - } - - template - ValueType const& LinearEquationSolver::getLowerBound() const { - return lowerBound.get(); - } - - template - ValueType const& LinearEquationSolver::getUpperBound() const { - return upperBound.get(); - } - - template - std::vector const& LinearEquationSolver::getLowerBounds() const { - return lowerBounds.get(); - } - - template - std::vector const& LinearEquationSolver::getUpperBounds() const { - return upperBounds.get(); - } - - template - void LinearEquationSolver::setLowerBounds(std::vector const& values) { - lowerBounds = values; - } - - template - void LinearEquationSolver::setUpperBounds(std::vector const& values) { - upperBounds = values; - } - - template - void LinearEquationSolver::setUpperBounds(std::vector&& values) { - upperBounds = std::move(values); - } - - template - void LinearEquationSolver::setBounds(std::vector const& lower, std::vector const& upper) { - setLowerBounds(lower); - setUpperBounds(upper); - } - - template - void LinearEquationSolver::createLowerBoundsVector(std::vector& lowerBoundsVector) const { - if (this->hasLowerBound(BoundType::Local)) { - lowerBoundsVector = this->getLowerBounds(); - } else { - STORM_LOG_THROW(this->hasLowerBound(BoundType::Global), storm::exceptions::UnmetRequirementException, "Cannot create lower bounds vector without lower bound."); - for (auto& e : lowerBoundsVector) { - e = this->getLowerBound(); - } - } - } - - template - void LinearEquationSolver::createUpperBoundsVector(std::unique_ptr>& upperBoundsVector) const { - if (!upperBoundsVector) { - if (this->hasUpperBound(BoundType::Local)) { - upperBoundsVector = std::make_unique>(this->getUpperBounds()); - } else { - upperBoundsVector = std::make_unique>(getMatrixRowCount(), this->getUpperBound()); - } - } else { - if (this->hasUpperBound(BoundType::Local)) { - for (auto& e : *upperBoundsVector) { - e = this->getUpperBound(); - } - } else { - auto upperBoundsIt = this->getUpperBounds().begin(); - for (auto& e : *upperBoundsVector) { - e = *upperBoundsIt; - ++upperBoundsIt; - } - } - } - } - template std::unique_ptr> LinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { std::unique_ptr> solver = this->create(); diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h index 780ae4ace..e2e5fdb11 100644 --- a/src/storm/solver/LinearEquationSolver.h +++ b/src/storm/solver/LinearEquationSolver.h @@ -150,98 +150,12 @@ namespace storm { */ virtual void clearCache() const; - enum class BoundType { - Global, - Local, - Any - }; - - /*! - * Retrieves whether this solver has a lower bound. - */ - bool hasLowerBound(BoundType const& type = BoundType::Any) const; - - /*! - * Retrieves whether this solver has an upper bound. - */ - bool hasUpperBound(BoundType const& type = BoundType::Any) const; - - /*! - * Sets a lower bound for the solution that can potentially be used by the solver. - */ - void setLowerBound(ValueType const& value); - - /*! - * Sets an upper bound for the solution that can potentially be used by the solver. - */ - void setUpperBound(ValueType const& value); - - /*! - * Sets bounds for the solution that can potentially be used by the solver. - */ - void setBounds(ValueType const& lower, ValueType const& upper); - - /*! - * Retrieves the lower bound (if there is any). - */ - ValueType const& getLowerBound() const; - - /*! - * Retrieves the upper bound (if there is any). - */ - ValueType const& getUpperBound() const; - - /*! - * Retrieves a vector containing the lower bounds (if there are any). - */ - std::vector const& getLowerBounds() const; - - /*! - * Retrieves a vector containing the upper bounds (if there are any). - */ - std::vector const& getUpperBounds() const; - - /*! - * Sets lower bounds for the solution that can potentially be used by the solver. - */ - void setLowerBounds(std::vector const& values); - - /*! - * Sets upper bounds for the solution that can potentially be used by the solver. - */ - void setUpperBounds(std::vector const& values); - - /*! - * Sets upper bounds for the solution that can potentially be used by the solver. - */ - void setUpperBounds(std::vector&& values); - - /*! - * Sets bounds for the solution that can potentially be used by the solver. - */ - void setBounds(std::vector const& lower, std::vector const& upper); - protected: virtual bool internalSolveEquations(std::vector& x, std::vector const& b) const = 0; - - void createUpperBoundsVector(std::unique_ptr>& upperBoundsVector) const; - void createLowerBoundsVector(std::vector& lowerBoundsVector) const; - + // auxiliary storage. If set, this vector has getMatrixRowCount() entries. mutable std::unique_ptr> cachedRowVector; - // A lower bound if one was set. - boost::optional lowerBound; - - // An upper bound if one was set. - boost::optional upperBound; - - // Lower bounds if they were set. - boost::optional> lowerBounds; - - // Lower bounds if they were set. - boost::optional> upperBounds; - private: /*! * Retrieves the row count of the matrix associated with this solver. diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index a61307047..06e690d0d 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -110,43 +110,7 @@ namespace storm { void MinMaxLinearEquationSolver::clearCache() const { // Intentionally left empty. } - - template - void MinMaxLinearEquationSolver::setLowerBound(ValueType const& value) { - lowerBound = value; - } - - template - void MinMaxLinearEquationSolver::setUpperBound(ValueType const& value) { - upperBound = value; - } - - template - void MinMaxLinearEquationSolver::setBounds(ValueType const& lower, ValueType const& upper) { - setLowerBound(lower); - setUpperBound(upper); - } - - template - bool MinMaxLinearEquationSolver::hasUpperBound() const { - return static_cast(upperBound); - } - - template - bool MinMaxLinearEquationSolver::hasLowerBound() const { - return static_cast(lowerBound); - } - - template - ValueType const& MinMaxLinearEquationSolver::getUpperBound() const { - return upperBound.get(); - } - - template - ValueType const& MinMaxLinearEquationSolver::getLowerBound() const { - return lowerBound.get(); - } - + template void MinMaxLinearEquationSolver::setInitialScheduler(std::vector&& choices) { initialScheduler = std::move(choices); diff --git a/src/storm/solver/MinMaxLinearEquationSolver.h b/src/storm/solver/MinMaxLinearEquationSolver.h index c05b2f13f..61742d326 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.h +++ b/src/storm/solver/MinMaxLinearEquationSolver.h @@ -135,41 +135,6 @@ namespace storm { * Clears the currently cached data that has been stored during previous calls of the solver. */ virtual void clearCache() const; - - /*! - * Sets a lower bound for the solution that can potentially used by the solver. - */ - void setLowerBound(ValueType const& value); - - /*! - * Sets an upper bound for the solution that can potentially used by the solver. - */ - void setUpperBound(ValueType const& value); - - /*! - * Sets bounds for the solution that can potentially used by the solver. - */ - void setBounds(ValueType const& lower, ValueType const& upper); - - /*! - * Retrieves whether the solver has an upper bound. - */ - bool hasUpperBound() const; - - /*! - * Retrieves whether the solver has a lower bound. - */ - bool hasLowerBound() const; - - /*! - * Retrieves the upper bound (if this solver has any). - */ - ValueType const& getUpperBound() const; - - /*! - * Retrieves the upper bound (if this solver has any). - */ - ValueType const& getLowerBound() const; /*! * Sets a valid initial scheduler that is required by some solvers (see requirements of solvers). @@ -205,7 +170,7 @@ namespace storm { protected: virtual bool internalSolveEquations(OptimizationDirection d, std::vector& x, std::vector const& b) const = 0; - + /// The optimization direction to use for calls to functions that do not provide it explicitly. Can also be unset. OptimizationDirectionSetting direction; @@ -215,12 +180,6 @@ namespace storm { /// The scheduler choices that induce the optimal values (if they could be successfully generated). mutable boost::optional> schedulerChoices; - // A lower bound if one was set. - boost::optional lowerBound; - - // An upper bound if one was set. - boost::optional upperBound; - // A scheduler that can be used by solvers that require a valid initial scheduler. boost::optional> initialScheduler; diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 2749c5219..b2db284c9 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -438,7 +438,7 @@ namespace storm { std::vector* lowerX = &x; this->createLowerBoundsVector(*lowerX); - this->createUpperBoundsVector(this->cachedRowVector); + this->createUpperBoundsVector(this->cachedRowVector, this->getMatrixRowCount()); std::vector* upperX = this->cachedRowVector.get(); bool useGaussSeidelMultiplication = this->getSettings().getPowerMethodMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel; diff --git a/src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h b/src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h index 922022da3..95491980a 100644 --- a/src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h +++ b/src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h @@ -2,6 +2,7 @@ #include #include +#include #include "storm/utility/macros.h"