fixed priority queue implementation and upper reward bound computation

src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp

@@ -6,6 +6,8 @@
 #include "storm/storage/BitVector.h"
 #include "storm/storage/StronglyConnectedComponentDecomposition.h"
 
+#include "storm/utility/macros.h"
+
 namespace storm {
     namespace modelchecker {
         namespace helper {
@@ -16,7 +18,9 @@ namespace storm {
             }
 
             template<typename ValueType>
-            std::vector<ValueType> BaierUpperRewardBoundsComputer<ValueType>::computeUpperBounds() {
+            ValueType BaierUpperRewardBoundsComputer<ValueType>::computeUpperBound() {
+                std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
+                
                 std::vector<uint64_t> stateToScc(transitionMatrix.getRowGroupCount());
                 {
                     // Start with an SCC decomposition of the system.
@@ -46,8 +50,11 @@ namespace storm {
                     ++index;
                 }
                 
+                // Vector that holds the result.
+                std::vector<ValueType> result(transitionMatrix.getRowGroupCount());
+                
                 // Process all states as long as there are remaining ones.
-                storm::storage::BitVector newStates(remainingStates.size());
+                std::vector<uint64_t> newStates;
                 while (!remainingStates.empty()) {
                     for (auto state : remainingStates) {
                         bool allChoicesValid = true;
@@ -55,45 +62,67 @@ namespace storm {
                             if (validChoices.get(row)) {
                                 continue;
                             }
+                            
+                            bool choiceValid = false;
                             for (auto const& entry : transitionMatrix.getRow(row)) {
                                 if (storm::utility::isZero(entry.getValue())) {
                                     continue;
                                 }
                                 
-                                if (remainingStates.get(entry.getColumn())) {
-                                    allChoicesValid = false;
+                                if (!remainingStates.get(entry.getColumn())) {
+                                    choiceValid = true;
                                     break;
                                 }
                             }
                             
-                            if (allChoicesValid) {
+                            if (choiceValid) {
                                 validChoices.set(row);
+                            } else {
+                                allChoicesValid = false;
                             }
                         }
                         
                         if (allChoicesValid) {
-                            newStates.set(state);
-                            remainingStates.set(state, false);
+                            newStates.push_back(state);
                         }
                     }
                     
                     // Compute d_t over the newly found states.
-                    ValueType d_t = storm::utility::one<ValueType>();
                     for (auto state : newStates) {
+                        result[state] = storm::utility::one<ValueType>();
                         for (auto row = transitionMatrix.getRowGroupIndices()[state], endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) {
-                            ValueType value = storm::utility::zero<ValueType>();
+                            ValueType rowValue = oneStepTargetProbabilities[row];
                             for (auto const& entry : transitionMatrix.getRow(row)) {
-                                if () {
-                                    
+                                if (!remainingStates.get(entry.getColumn())) {
+                                    rowValue += entry.getValue() * (stateToScc[state] == stateToScc[entry.getColumn()] ? result[entry.getColumn()] : storm::utility::one<ValueType>());
                                 }
                             }
-                            d_t = std::min();
+                            STORM_LOG_ASSERT(rowValue > storm::utility::zero<ValueType>(), "Expected entry with value greater 0.");
+                            result[state] = std::min(result[state], rowValue);
                         }
                     }
+                    
+                    remainingStates.set(newStates.begin(), newStates.end(), false);
                     newStates.clear();
                 }
                 
+                for (uint64_t state = 0; state < result.size(); ++state) {
+                    ValueType maxReward = storm::utility::zero<ValueType>();
+                    for (auto row = transitionMatrix.getRowGroupIndices()[state], endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) {
+                        maxReward = std::max(maxReward, rewards[row]);
+                    }
+                    result[state] = storm::utility::one<ValueType>() / result[state] * maxReward;
+                }
+                
+                ValueType upperBound = storm::utility::zero<ValueType>();
+                for (auto const& e : result) {
+                    upperBound += e;
+                }
                 
+                STORM_LOG_TRACE("Baier algorithm for reward bound computation (variant 2) computed bound " << upperBound << ".");
+                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 upperBound;
             }
             
             template class BaierUpperRewardBoundsComputer<double>;

src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h

@@ -25,9 +25,9 @@ namespace storm {
                 BaierUpperRewardBoundsComputer(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> const& oneStepTargetProbabilities);
                 
                 /*!
-                 * Computes upper bounds on the expected rewards.
+                 * Computes an upper bound on the expected rewards.
                  */
-                std::vector<ValueType> computeUpperBounds();
+                ValueType computeUpperBound();
                 
             private:
                 storm::storage::SparseMatrix<ValueType> const& transitionMatrix;

src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp

@@ -23,6 +23,7 @@ namespace storm {
             
             template<typename ValueType>
             std::vector<ValueType> DsMpiDtmcUpperRewardBoundsComputer<ValueType>::computeUpperBounds() {
+                std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
                 sweep();
                 ValueType lambda = computeLambda();
                 STORM_LOG_TRACE("DS-MPI computed lambda as " << lambda << ".");
@@ -45,6 +46,8 @@ namespace storm {
                 }
                 STORM_LOG_TRACE("DS-MPI computed " << nonZeroCount << " non-zero upper bounds and a maximal bound of " << max << ".");
 #endif
+                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 result;
             }
             
@@ -63,16 +66,15 @@ namespace storm {
                 uint64_t state = this->getStateForChoice(choice);
                 
                 // Check whether condition (I) or (II) applies.
-                ValueType sum = storm::utility::zero<ValueType>();
+                ValueType probSum = originalOneStepTargetProbabilities[choice];
                 for (auto const& e : transitionMatrix.getRow(choice)) {
-                    sum += e.getValue() * p[e.getColumn()];
+                    probSum += 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 << ".");
+                if (p[state] < probSum) {
+                    STORM_LOG_TRACE("Condition (I) does apply for state " << state << " as " << p[state] << " < " << probSum << ".");
                     // Condition (I) applies.
-                    localLambda = sum - p[state];
+                    localLambda = probSum - p[state];
                     ValueType nominator = originalRewards[choice];
                     for (auto const& e : transitionMatrix.getRow(choice)) {
                         nominator += e.getValue() * w[e.getColumn()];
@@ -80,17 +82,18 @@ namespace storm {
                     nominator -= w[state];
                     localLambda = nominator / localLambda;
                 } else {
-                    STORM_LOG_TRACE("Condition (I) does not apply for state " << state << " as " << p[state] << " < " << sum << ".");
+                    STORM_LOG_TRACE("Condition (I) does not apply for state " << state << std::setprecision(30) << " as " << probSum << " <= " << p[state] << ".");
                     // 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];
+                    ValueType rewardSum = originalRewards[choice];
                     for (auto const& e : transitionMatrix.getRow(choice)) {
-                        sum += e.getValue() * w[e.getColumn()];
+                        rewardSum += e.getValue() * w[e.getColumn()];
                     }
-                    STORM_LOG_WARN_COND(w[state] >= sum || storm::utility::ConstantsComparator<ValueType>().isEqual(w[state], sum), "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << sum << ".");
+                    STORM_LOG_WARN_COND(w[state] >= rewardSum || storm::utility::ConstantsComparator<ValueType>().isEqual(w[state], rewardSum), "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << rewardSum << ".");
+                    STORM_LOG_WARN_COND(storm::utility::ConstantsComparator<ValueType>().isEqual(probSum, p[state]), "Expected condition (II) to hold in state " << state << ", but " << probSum << " != " << p[state] << ".");
 #endif
                 }
 
@@ -237,7 +240,7 @@ namespace storm {
                 while (!queue.empty()) {
                     // Get first entry in queue.
                     storm::storage::sparse::state_type currentState = queue.popTop();
-                    
+
                     // Mark state as visited.
                     visited.set(currentState);
                     
@@ -245,7 +248,7 @@ namespace storm {
                     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)) {

src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp

@@ -224,11 +224,8 @@ 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();
                 DsMpiDtmcUpperRewardBoundsComputer<ValueType> dsmpi(transitionMatrix, rewards, oneStepTargetProbabilities);
                 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;
             }
             

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -5,6 +5,7 @@
 #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "storm/modelchecker/hints/ExplicitModelCheckerHint.h"
 #include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h"
+#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h"
 
 #include "storm/models/sparse/StandardRewardModel.h"
 
@@ -1020,7 +1021,8 @@ namespace storm {
                     DsMpiMdpUpperRewardBoundsComputer<ValueType> dsmpi(submatrix, choiceRewards, oneStepTargetProbabilities);
                     hintInformation.upperResultBounds = dsmpi.computeUpperBounds();
                 } else {
-                    STORM_LOG_ASSERT(false, "Not yet implemented.");
+                    BaierUpperRewardBoundsComputer<ValueType> baier(submatrix, choiceRewards, oneStepTargetProbabilities);
+                    hintInformation.upperResultBound = baier.computeUpperBound();
                 }
             }
             
@@ -1079,7 +1081,7 @@ namespace storm {
                             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, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b);
+                            computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr);
                         }
                         
                         // If we need to compute upper bounds, do so now.

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -410,7 +410,7 @@ namespace storm {
                 }
                 
                 // Determine whether the method converged.
-                if (storm::utility::vector::equalModuloPrecision<ValueType>(*lowerX, *upperX, storm::utility::convertNumber<ValueType>(2.0) * this->getSettings().getPrecision(), false)) {
+                if (storm::utility::vector::equalModuloPrecision<ValueType>(*lowerX, *upperX, storm::utility::convertNumber<ValueType>(2.0) * this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) {
                     status = Status::Converged;
                 }
                 

src/storm/storage/BitVector.cpp

@@ -184,9 +184,9 @@ namespace storm {
         }
 
         template<typename InputIterator>
-        void BitVector::set(InputIterator begin, InputIterator end) {
+        void BitVector::set(InputIterator begin, InputIterator end, bool value) {
             for (InputIterator it = begin; it != end; ++it) {
-                this->set(*it);
+                this->set(*it, value);
             }
         }
 
@@ -1020,10 +1020,10 @@ namespace storm {
         template BitVector::BitVector(uint_fast64_t length, std::vector<uint_fast64_t>::const_iterator begin, std::vector<uint_fast64_t>::const_iterator end);
         template BitVector::BitVector(uint_fast64_t length, boost::container::flat_set<uint_fast64_t>::iterator begin, boost::container::flat_set<uint_fast64_t>::iterator end);
         template BitVector::BitVector(uint_fast64_t length, boost::container::flat_set<uint_fast64_t>::const_iterator begin, boost::container::flat_set<uint_fast64_t>::const_iterator end);
-        template void BitVector::set(std::vector<uint_fast64_t>::iterator begin, std::vector<uint_fast64_t>::iterator end);
-        template void BitVector::set(std::vector<uint_fast64_t>::const_iterator begin, std::vector<uint_fast64_t>::const_iterator end);
-        template void BitVector::set(boost::container::flat_set<uint_fast64_t>::iterator begin, boost::container::flat_set<uint_fast64_t>::iterator end);
-        template void BitVector::set(boost::container::flat_set<uint_fast64_t>::const_iterator begin, boost::container::flat_set<uint_fast64_t>::const_iterator end);
+        template void BitVector::set(std::vector<uint_fast64_t>::iterator begin, std::vector<uint_fast64_t>::iterator end, bool value);
+        template void BitVector::set(std::vector<uint_fast64_t>::const_iterator begin, std::vector<uint_fast64_t>::const_iterator end, bool value);
+        template void BitVector::set(boost::container::flat_set<uint_fast64_t>::iterator begin, boost::container::flat_set<uint_fast64_t>::iterator end, bool value);
+        template void BitVector::set(boost::container::flat_set<uint_fast64_t>::const_iterator begin, boost::container::flat_set<uint_fast64_t>::const_iterator end, bool value);
     }
 }
 

src/storm/storage/BitVector.h

@@ -205,7 +205,7 @@ namespace storm {
              * @param last The element past the last element of the iterator range.
              */
             template<typename InputIterator>
-            void set(InputIterator first, InputIterator last);
+            void set(InputIterator first, InputIterator last, bool value = true);
             
             /*!
              * Retrieves the truth value of the bit at the given index. Note: this does not check whether the given

src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h

@@ -4,6 +4,8 @@
 #include <vector>
 #include <numeric>
 
+#include "storm-config.h"
+
 #include "storm/utility/macros.h"
 
 namespace storm {
@@ -24,10 +26,8 @@ namespace storm {
         public:
             explicit ConsecutiveUint64DynamicPriorityQueue(uint64_t numberOfIntegers, Compare const& compare) : container(numberOfIntegers), compare(compare), positions(numberOfIntegers) {
                 std::iota(container.begin(), container.end(), 0);
-            }
-            
-            void fix() {
                 std::make_heap(container.begin(), container.end(), compare);
+                updatePositions();
             }
             
             void increase(uint64_t element) {
@@ -38,12 +38,14 @@ namespace storm {
                 
                 uint64_t parentPosition = (position - 1) / 2;
                 while (position > 0 && compare(container[parentPosition], container[position])) {
-                    std::swap(container[parentPosition], container[position]);
                     std::swap(positions[container[parentPosition]], positions[container[position]]);
-                    
+                    std::swap(container[parentPosition], container[position]);
+
                     position = parentPosition;
                     parentPosition = (position - 1) / 2;
                 }
+                
+                STORM_LOG_ASSERT(std::is_heap(container.begin(), container.end(), compare), "Heap structure lost.");
             }
             
             bool contains(uint64_t element) const {
@@ -68,8 +70,47 @@ namespace storm {
             }
             
             void pop() {
-                std::pop_heap(container.begin(), container.end(), compare);
-                container.pop_back();
+                if (container.size() > 1) {
+                    // Swap max element to back.
+                    std::swap(positions[container.front()], positions[container.back()]);
+                    std::swap(container.front(), container.back());
+                    container.pop_back();
+                    
+                    // Sift down the element from the top.
+                    uint64_t positionToSift = 0;
+                    uint64_t child = 2 * positionToSift + 1;
+                    
+                    while (child < container.size()) {
+                        if (child + 1 < container.size()) {
+                            // Figure out larger child.
+                            child = compare(container[child], container[child + 1]) ? child + 1 : child;
+                            
+                            // Check if we need to sift down.
+                            if (compare(container[positionToSift], container[child])) {
+                                std::swap(positions[container[positionToSift]], positions[container[child]]);
+                                std::swap(container[positionToSift], container[child]);
+                                
+                                positionToSift = child;
+                                child = 2 * positionToSift + 1;
+                            } else {
+                                break;
+                            }
+                        } else if (compare(container[positionToSift], container[child])) {
+                            std::swap(positions[container[positionToSift]], positions[container[child]]);
+                            std::swap(container[positionToSift], container[child]);
+                            
+                            positionToSift = child;
+                            child = 2 * positionToSift + 1;
+                        } else {
+                            break;
+                        }
+                    }
+                    
+                } else {
+                    container.pop_back();
+                }
+                
+                STORM_LOG_ASSERT(std::is_heap(container.begin(), container.end(), compare), "Heap structure lost.");
             }
             
             T popTop() {
@@ -77,6 +118,26 @@ namespace storm {
                 pop();
                 return item;
             }
+            
+        private:
+            bool checkPositions() const {
+                uint64_t position = 0;
+                for (auto const& e : container) {
+                    if (positions[e] != position) {
+                        return false;
+                    }
+                    ++position;
+                }
+                return true;
+            }
+            
+            void updatePositions() {
+                uint64_t position = 0;
+                for (auto const& e : container) {
+                    positions[e] = position;
+                    ++position;
+                }
+            }
         };
     }
 }