diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 90ed634db..0635c54cd 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -7,6 +7,7 @@ #include "storm/utility/graph.h" #include "storm/storage/StronglyConnectedComponentDecomposition.h" +#include "storm/storage/DynamicPriorityQueue.h" #include "storm/solver/LinearEquationSolver.h" @@ -17,6 +18,7 @@ #include "storm/exceptions/InvalidPropertyException.h" #include "storm/exceptions/IllegalArgumentException.h" #include "storm/exceptions/UncheckedRequirementException.h" +#include "storm/exceptions/NotSupportedException.h" namespace storm { namespace modelchecker { @@ -217,13 +219,14 @@ namespace storm { template class DsMpi { public: - DsMpi(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) : transitionMatrix(transitionMatrix), originalRewards(rewards), backwardTransitions(transitionMatrix.transpose()), p(transitionMatrix.getRowCount()), w(transitionMatrix.getRowCount()), rewards(rewards), targetProbabilities(oneStepTargetProbabilities) { + 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()); @@ -231,18 +234,21 @@ namespace storm { for (storm::storage::sparse::state_type state = 0; state < result.size(); ++state) { result[state] = w[state] + (one - p[state]) * lambda; } + return result; } private: ValueType computeLambda() { - ValueType lambda = storm::utility::convertNumber(0.0); + 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]; @@ -256,38 +262,59 @@ namespace storm { } 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, "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << sum << "."); +#endif } } return lambda; } - void sweep() { - // Create a priority queue that allows for easy retrieval of the currently best state. - auto cmp = [this](storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { - ValueType pa = p[a]; - ValueType pb = p[b]; - if (pa > pb) { + 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 w[a] < w[b]; + return dsmpi.rewards[a] > dsmpi.rewards[b]; } return false; - }; - std::set queue; - + } + + private: + DsMpi const& dsmpi; + }; + + void sweep() { + // Create a priority queue that allows for easy retrieval of the currently best state. + storm::storage::DynamicPriorityQueue, PriorityLess> queue(PriorityLess(*this)); storm::storage::BitVector visited(p.size()); + storm::storage::BitVector inQueue(p.size()); for (storm::storage::sparse::state_type state = 0; state < targetProbabilities.size(); ++state) { if (!storm::utility::isZero(targetProbabilities[state])) { - queue.insert(state); + queue.push(state); + inQueue.set(state); } } + queue.fix(); while (!queue.empty()) { // Get first entry in queue. - storm::storage::sparse::state_type currentState = *queue.begin(); - queue.erase(queue.begin()); + storm::storage::sparse::state_type currentState = queue.popTop(); // Mark state as visited. visited.set(currentState); @@ -301,18 +328,17 @@ namespace storm { continue; } - // Delete element from the priority queue if it was in there. - auto it = queue.find(e.getColumn()); - if (it != queue.end()) { - queue.erase(it); - } - // Update reward/probability values. rewards[e.getColumn()] += e.getValue() * w[currentState]; targetProbabilities[e.getColumn()] += e.getValue() * p[currentState]; - // (Re-)insert the state with the new rewards/target probabilities. - queue.insert(e.getColumn()); + // Either insert element or simply fix the queue. + if (!inQueue.get(e.getColumn())) { + queue.push(e.getColumn()); + inQueue.set(e.getColumn()); + } else { + queue.fix(); + } } } } @@ -320,6 +346,7 @@ namespace storm { // 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; @@ -334,8 +361,17 @@ namespace storm { // 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); - return dsmpi.computeUpperBounds(); + 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."); + return bounds; + } + + template<> + std::vector computeUpperRewardBounds(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& rewards, std::vector const& oneStepTargetProbabilities) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing upper reward bounds is not supported for rational functions."); } template diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index 88e06293c..22c75794b 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -15,6 +15,7 @@ #include "storm/utility/macros.h" #include "storm/exceptions/NotSupportedException.h" +#include "storm/exceptions/UnmetRequirementException.h" namespace storm { namespace solver { @@ -212,6 +213,18 @@ namespace storm { 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) { diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h index dc80bd688..780ae4ace 100644 --- a/src/storm/solver/LinearEquationSolver.h +++ b/src/storm/solver/LinearEquationSolver.h @@ -225,6 +225,7 @@ namespace storm { 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; diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index c4cdab06c..393e44be1 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -382,6 +382,7 @@ namespace storm { template bool NativeLinearEquationSolver::solveEquationsPower(std::vector& x, std::vector const& b) const { + STORM_LOG_THROW(this->hasLowerBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Power)"); if (!this->cachedRowVector) { @@ -389,6 +390,7 @@ namespace storm { } std::vector* currentX = &x; + this->createLowerBoundsVector(*currentX); std::vector* nextX = this->cachedRowVector.get(); bool useGaussSeidelMultiplication = this->getSettings().getPowerMethodMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel; @@ -430,10 +432,12 @@ namespace storm { template bool NativeLinearEquationSolver::solveEquationsSoundPower(std::vector& x, std::vector const& b) const { + STORM_LOG_THROW(this->hasLowerBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); STORM_LOG_THROW(this->hasUpperBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (SoundPower)"); std::vector* lowerX = &x; + this->createLowerBoundsVector(*lowerX); this->createUpperBoundsVector(this->cachedRowVector); std::vector* upperX = this->cachedRowVector.get(); @@ -446,24 +450,65 @@ namespace storm { bool converged = false; uint64_t iterations = 0; + bool doConvergenceCheck = false; + ValueType upperDiff; + ValueType lowerDiff; while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations()) { - if (useGaussSeidelMultiplication) { - this->multiplier.multAddGaussSeidelBackward(*this->A, *lowerX, &b); - this->multiplier.multAddGaussSeidelBackward(*this->A, *upperX, &b); + // In every hundredth iteration, we improve both bounds. + if (iterations % 100 == 0) { + if (useGaussSeidelMultiplication) { + lowerDiff = (*lowerX)[0]; + this->multiplier.multAddGaussSeidelBackward(*this->A, *lowerX, &b); + lowerDiff = (*lowerX)[0] - lowerDiff; + upperDiff = (*upperX)[0]; + this->multiplier.multAddGaussSeidelBackward(*this->A, *upperX, &b); + upperDiff = upperDiff - (*upperX)[0]; + } else { + this->multiplier.multAdd(*this->A, *lowerX, &b, *tmp); + lowerDiff = (*tmp)[0] - (*lowerX)[0]; + std::swap(tmp, lowerX); + this->multiplier.multAdd(*this->A, *upperX, &b, *tmp); + upperDiff = (*upperX)[0] - (*tmp)[0]; + std::swap(tmp, upperX); + } } else { - this->multiplier.multAdd(*this->A, *lowerX, &b, *tmp); - std::swap(tmp, lowerX); - this->multiplier.multAdd(*this->A, *upperX, &b, *tmp); - std::swap(tmp, upperX); + // In the following iterations, we improve the bound with the greatest difference. + if (useGaussSeidelMultiplication) { + if (lowerDiff >= upperDiff) { + lowerDiff = (*lowerX)[0]; + this->multiplier.multAddGaussSeidelBackward(*this->A, *lowerX, &b); + lowerDiff = (*lowerX)[0] - lowerDiff; + } else { + upperDiff = (*upperX)[0]; + this->multiplier.multAddGaussSeidelBackward(*this->A, *upperX, &b); + upperDiff = upperDiff - (*upperX)[0]; + } + } else { + if (lowerDiff >= upperDiff) { + this->multiplier.multAdd(*this->A, *lowerX, &b, *tmp); + lowerDiff = (*tmp)[0] - (*lowerX)[0]; + std::swap(tmp, lowerX); + } else { + this->multiplier.multAdd(*this->A, *upperX, &b, *tmp); + upperDiff = (*upperX)[0] - (*tmp)[0]; + std::swap(tmp, upperX); + } + } } + STORM_LOG_ASSERT(lowerDiff >= storm::utility::zero(), "Expected non-negative lower diff."); + STORM_LOG_ASSERT(upperDiff >= storm::utility::zero(), "Expected non-negative upper diff."); + STORM_LOG_TRACE("Lower difference: " << lowerDiff << ", upper difference: " << upperDiff << "."); - // Now check if the process already converged within our precision. Note that we double the target - // precision here. Doing so, we need to take the means of the lower and upper values later to guarantee - // the original precision. - converged = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, storm::utility::convertNumber(2.0) * static_cast(this->getSettings().getPrecision()), false); + if (doConvergenceCheck) { + // Now check if the process already converged within our precision. Note that we double the target + // precision here. Doing so, we need to take the means of the lower and upper values later to guarantee + // the original precision. + converged = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, storm::utility::convertNumber(2.0) * static_cast(this->getSettings().getPrecision()), false); + } // Set up next iteration. ++iterations; + doConvergenceCheck = !doConvergenceCheck; } // We take the means of the lower and upper bound so we guarantee the desired precision. @@ -582,6 +627,24 @@ namespace storm { } } + template + LinearEquationSolverRequirements NativeLinearEquationSolver::getRequirements() const { + LinearEquationSolverRequirements requirements; + if (this->getSettings().getForceSoundness()) { + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::Power) { + requirements.requireLowerBounds(); + requirements.requireUpperBounds(); + } else { + STORM_LOG_WARN("Forcing soundness, but selecting a method other than the power iteration is not supported."); + } + } else { + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::Power) { + requirements.requireLowerBounds(); + } + } + return requirements; + } + template void NativeLinearEquationSolver::clearCache() const { jacobiDecomposition.reset(); diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index c252b86b3..d0c2ae410 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -68,6 +68,7 @@ namespace storm { NativeLinearEquationSolverSettings const& getSettings() const; virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override; + virtual LinearEquationSolverRequirements getRequirements() const override; virtual void clearCache() const override; diff --git a/src/storm/storage/DynamicPriorityQueue.h b/src/storm/storage/DynamicPriorityQueue.h index 9df114fbf..d54f7c977 100644 --- a/src/storm/storage/DynamicPriorityQueue.h +++ b/src/storm/storage/DynamicPriorityQueue.h @@ -9,7 +9,6 @@ namespace storm { template, typename Compare = std::less> class DynamicPriorityQueue { - private: Container container; Compare compare; diff --git a/src/storm/storage/VectorDynamicPriorityQueue.h b/src/storm/storage/VectorDynamicPriorityQueue.h new file mode 100644 index 000000000..46931d3a3 --- /dev/null +++ b/src/storm/storage/VectorDynamicPriorityQueue.h @@ -0,0 +1,99 @@ +#ifndef STORM_STORAGE_DYNAMICPRIORITYQUEUE_H_ +#define STORM_STORAGE_DYNAMICPRIORITYQUEUE_H_ + +#include +#include + +#include "storm/utility/macros.h" + +namespace storm { + namespace storage { + + template> + class VectorDynamicPriorityQueue { + public: + typedef uint64_t T; + typedef std::vector Container; + + private: + Container container; + Compare compare; + + std::vector positions; + uint64_t upperBound; + + uint64_t numberOfSortedEntriesAtBack; + + uint64_t const NUMBER_OF_ENTRIES_TO_SORT = 100; + + public: + explicit DynamicPriorityQueue(Compare const& compare, uint64_t upperBound) : container(), compare(compare), positions(upperBound) { + // Intentionally left empty + } + + explicit DynamicPriorityQueue(Container&& container, Compare const& compare) : container(std::move(container)), compare(compare), positions(this->container.size()) { + sortAndUpdatePositions(container.begin(), container.end()); + } + + void fix() { + sortAndUpdatePositions(container.begin(), container.end()); + } + + bool empty() const { + return container.empty(); + } + + std::size_t size() const { + return container.size(); + } + + const T& top() const { + return container.front(); + } + + template + void push(TemplateType&& item) { + if (this->empty() || container.back() < item) { + container.emplace_back(std::forward(item)); + } else { + + } + } + + void pop() { + container.pop_back(); + --numberOfSortedEntriesAtBack; + if (numberOfSortedEntriesAtBack == 0) { + if (container.size() > NUMBER_OF_ENTRIES_TO_SORT) { + sortAndUpdatePositions(container.end() - NUMBER_OF_ENTRIES_TO_SORT, container.end()); + numberOfSortedEntriesAtBack = NUMBER_OF_ENTRIES_TO_SORT; + } else { + sortAndUpdatePositions(container.begin(), container.end()); + numberOfSortedEntriesAtBack = container.size(); + } + } + } + + T popTop() { + T item = top(); + pop(); + return item; + } + + private: + void sortAndUpdatePositions(Container::const_iterator start, Container::const_iterator end) { + std::sort(start, end); + updatePositions(start, end); + } + + void updatePositions(Container::const_iterator start, Container::const_iterator end) { + for (; start != end; ++start) { + position = std::distance(container.begin(), start); + positions[container[position]] = position; + } + } + }; + } +} + +#endif // STORM_STORAGE_DYNAMICPRIORITYQUEUE_H_