intermediate commit

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<typename ValueType>
             class DsMpi {
             public:
-                DsMpi(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> const& oneStepTargetProbabilities) : transitionMatrix(transitionMatrix), originalRewards(rewards), backwardTransitions(transitionMatrix.transpose()), p(transitionMatrix.getRowCount()), w(transitionMatrix.getRowCount()), rewards(rewards), targetProbabilities(oneStepTargetProbabilities) {
+                DsMpi(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> 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<ValueType> computeUpperBounds() {
                     sweep();
                     ValueType lambda = computeLambda();
+                    STORM_LOG_TRACE("DS-MPI computed lambda as " << lambda << ".");
                     
                     // Finally compute the upper bounds for the states.
                     std::vector<ValueType> 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<ValueType>(0.0);
+                    ValueType lambda = storm::utility::zero<ValueType>();
                     for (storm::storage::sparse::state_type state = 0; state < targetProbabilities.size(); ++state) {
                         // Check whether condition (I) or (II) applies.
                         ValueType sum = storm::utility::zero<ValueType>();
                         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<storm::storage::sparse::state_type, decltype(cmp)> queue;
+                    }
                     
+                private:
+                    DsMpi const& dsmpi;
+                };
+                
+                void sweep() {
+                    // Create a priority queue that allows for easy retrieval of the currently best state.
+                    storm::storage::DynamicPriorityQueue<storm::storage::sparse::state_type, std::vector<storm::storage::sparse::state_type>, 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<ValueType> const& transitionMatrix;
                 std::vector<ValueType> const& originalRewards;
+                std::vector<ValueType> const& originalOneStepTargetProbabilities;
                 
                 // Derived from input data.
                 storm::storage::SparseMatrix<ValueType> backwardTransitions;
@@ -334,8 +361,17 @@ namespace storm {
             // This function computes an upper bound on the reachability rewards (see Baier et al, CAV'17).
             template<typename ValueType>
             std::vector<ValueType> computeUpperRewardBounds(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> const& oneStepTargetProbabilities) {
+                std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
                 DsMpi<ValueType> dsmpi(transitionMatrix, rewards, oneStepTargetProbabilities);
-                return dsmpi.computeUpperBounds();
+                std::vector<ValueType> 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<std::chrono::milliseconds>(end - start).count() << "ms.");
+                return bounds;
+            }
+            
+            template<>
+            std::vector<storm::RationalFunction> computeUpperRewardBounds(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, std::vector<storm::RationalFunction> const& rewards, std::vector<storm::RationalFunction> const& oneStepTargetProbabilities) {
+                STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing upper reward bounds is not supported for rational functions.");
             }
             
             template<typename ValueType, typename RewardModelType>

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<typename ValueType>
+        void LinearEquationSolver<ValueType>::createLowerBoundsVector(std::vector<ValueType>& 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<typename ValueType>
         void LinearEquationSolver<ValueType>::createUpperBoundsVector(std::unique_ptr<std::vector<ValueType>>& upperBoundsVector) const {
             if (!upperBoundsVector) {

src/storm/solver/LinearEquationSolver.h

@@ -225,6 +225,7 @@ namespace storm {
             virtual bool internalSolveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const = 0;
             
             void createUpperBoundsVector(std::unique_ptr<std::vector<ValueType>>& upperBoundsVector) const;
+            void createLowerBoundsVector(std::vector<ValueType>& lowerBoundsVector) const;
             
             // auxiliary storage. If set, this vector has getMatrixRowCount() entries.
             mutable std::unique_ptr<std::vector<ValueType>> cachedRowVector;

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -382,6 +382,7 @@ namespace storm {
         
         template<typename ValueType>
         bool NativeLinearEquationSolver<ValueType>::solveEquationsPower(std::vector<ValueType>& x, std::vector<ValueType> 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<ValueType>* currentX = &x;
+            this->createLowerBoundsVector(*currentX);
             std::vector<ValueType>* nextX = this->cachedRowVector.get();
 
             bool useGaussSeidelMultiplication = this->getSettings().getPowerMethodMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;
@@ -430,10 +432,12 @@ namespace storm {
         
         template<typename ValueType>
         bool NativeLinearEquationSolver<ValueType>::solveEquationsSoundPower(std::vector<ValueType>& x, std::vector<ValueType> 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<ValueType>* lowerX = &x;
+            this->createLowerBoundsVector(*lowerX);
             this->createUpperBoundsVector(this->cachedRowVector);
             std::vector<ValueType>* 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<ValueType>(), "Expected non-negative lower diff.");
+                STORM_LOG_ASSERT(upperDiff >= storm::utility::zero<ValueType>(), "Expected non-negative upper diff.");
+                STORM_LOG_TRACE("Lower difference: " << lowerDiff << ", upper difference: " << upperDiff << ".");
+                
+                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<ValueType>(*lowerX, *upperX, storm::utility::convertNumber<ValueType>(2.0) * static_cast<ValueType>(this->getSettings().getPrecision()), false);
                 }
-                
-                // 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<ValueType>(*lowerX, *upperX, storm::utility::convertNumber<ValueType>(2.0) * static_cast<ValueType>(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<typename ValueType>
+        LinearEquationSolverRequirements NativeLinearEquationSolver<ValueType>::getRequirements() const {
+            LinearEquationSolverRequirements requirements;
+            if (this->getSettings().getForceSoundness()) {
+                if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::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<ValueType>::SolutionMethod::Power) {
+                    requirements.requireLowerBounds();
+                }
+            }
+            return requirements;
+        }
+        
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::clearCache() const {
             jacobiDecomposition.reset();

src/storm/solver/NativeLinearEquationSolver.h

@@ -68,6 +68,7 @@ namespace storm {
             NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
 
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
+            virtual LinearEquationSolverRequirements getRequirements() const override;
 
             virtual void clearCache() const override;
 

src/storm/storage/DynamicPriorityQueue.h

@@ -9,7 +9,6 @@ namespace storm {
 
         template<typename T, typename Container = std::vector<T>, typename Compare = std::less<T>>
         class DynamicPriorityQueue {
-
         private:
             Container container;
             Compare compare;

src/storm/storage/VectorDynamicPriorityQueue.h

@@ -0,0 +1,99 @@
+#ifndef STORM_STORAGE_DYNAMICPRIORITYQUEUE_H_
+#define STORM_STORAGE_DYNAMICPRIORITYQUEUE_H_
+
+#include <algorithm>
+#include <vector>
+
+#include "storm/utility/macros.h"
+
+namespace storm {
+    namespace storage {
+        
+        template<typename Compare = std::less<uint64_t>>
+        class VectorDynamicPriorityQueue {
+        public:
+            typedef uint64_t T;
+            typedef std::vector<T> Container;
+            
+        private:
+            Container container;
+            Compare compare;
+            
+            std::vector<uint64_t> 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<typename TemplateType>
+            void push(TemplateType&& item) {
+                if (this->empty() || container.back() < item) {
+                    container.emplace_back(std::forward<TemplateType>(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_