upper bounds computation for reachability rewards in sparse MDPs

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<typename ValueType>
+            DsMpiDtmcUpperRewardBoundsComputer<ValueType>::DsMpiDtmcUpperRewardBoundsComputer(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.getRowGroupCount()), w(transitionMatrix.getRowGroupCount()), rewards(rewards), targetProbabilities(oneStepTargetProbabilities) {
+                // Intentionally left empty.
+            }
+            
+            template<typename ValueType>
+            std::vector<ValueType> DsMpiDtmcUpperRewardBoundsComputer<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.getRowGroupCount());
+                auto one = storm::utility::one<ValueType>();
+                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<ValueType>();
+                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<typename ValueType>
+            ValueType DsMpiDtmcUpperRewardBoundsComputer<ValueType>::computeLambda() const {
+                ValueType lambda = storm::utility::zero<ValueType>();
+                for (storm::storage::sparse::state_type state = 0; state < transitionMatrix.getRowGroupCount(); ++state) {
+                    lambda = std::max(lambda, computeLambdaForChoice(state));
+                }
+                return lambda;
+            }
+            
+            template<typename ValueType>
+            ValueType DsMpiDtmcUpperRewardBoundsComputer<ValueType>::computeLambdaForChoice(uint64_t choice) const {
+                ValueType localLambda = storm::utility::zero<ValueType>();
+                uint64_t state = this->getStateForChoice(choice);
+                
+                // Check whether condition (I) or (II) applies.
+                ValueType sum = storm::utility::zero<ValueType>();
+                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<ValueType>().isEqual(w[state], sum), "Expected condition (II) to hold in state " << state << ", but " << w[state] << " < " << sum << ".");
+#endif
+                }
+
+                return localLambda;
+            }
+            
+            template<typename ValueType>
+            uint64_t DsMpiDtmcUpperRewardBoundsComputer<ValueType>::getStateForChoice(uint64_t choice) const {
+                return choice;
+            }
+            
+            template<typename ValueType>
+            class DsMpiDtmcPriorityLess {
+            public:
+                DsMpiDtmcPriorityLess(DsMpiDtmcUpperRewardBoundsComputer<ValueType> 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<ValueType> const& dsmpi;
+            };
+            
+            template<typename ValueType>
+            void DsMpiDtmcUpperRewardBoundsComputer<ValueType>::sweep() {
+                // Create a priority queue that allows for easy retrieval of the currently best state.
+                storm::storage::ConsecutiveUint64DynamicPriorityQueue<DsMpiDtmcPriorityLess<ValueType>> queue(transitionMatrix.getRowCount(), DsMpiDtmcPriorityLess<ValueType>(*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<typename ValueType>
+            DsMpiMdpUpperRewardBoundsComputer<ValueType>::DsMpiMdpUpperRewardBoundsComputer(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> const& oneStepTargetProbabilities) : DsMpiDtmcUpperRewardBoundsComputer<ValueType>(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<ValueType> minReward;
+                    ValueType maxProb = storm::utility::zero<ValueType>();
+                    
+                    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<typename ValueType>
+            ValueType DsMpiMdpUpperRewardBoundsComputer<ValueType>::computeLambda() const {
+                ValueType lambda = storm::utility::zero<ValueType>();
+                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<typename ValueType>
+            uint64_t DsMpiMdpUpperRewardBoundsComputer<ValueType>::getStateForChoice(uint64_t choice) const {
+                return choiceToState[choice];
+            }
+            
+            template<typename ValueType>
+            class DsMpiMdpPriorityLess {
+            public:
+                DsMpiMdpPriorityLess(DsMpiMdpUpperRewardBoundsComputer<ValueType> 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<ValueType> const& dsmpi;
+            };
+            
+            template<typename ValueType>
+            void DsMpiMdpUpperRewardBoundsComputer<ValueType>::sweep() {
+                // Create a priority queue that allows for easy retrieval of the currently best state.
+                storm::storage::ConsecutiveUint64DynamicPriorityQueue<DsMpiMdpPriorityLess<ValueType>> queue(this->transitionMatrix.getRowGroupCount(), DsMpiMdpPriorityLess<ValueType>(*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<typename ValueType>
+            uint64_t DsMpiMdpUpperRewardBoundsComputer<ValueType>::getChoiceInState(uint64_t state) const {
+                return policy[state];
+            }
+
+            template<typename ValueType>
+            void DsMpiMdpUpperRewardBoundsComputer<ValueType>::setChoiceInState(uint64_t state, uint64_t choice) {
+                policy[state] = choice;
+            }
+
+            template class DsMpiDtmcUpperRewardBoundsComputer<double>;
+            template class DsMpiMdpUpperRewardBoundsComputer<double>;
+
+#ifdef STORM_HAVE_CARL
+            template class DsMpiDtmcUpperRewardBoundsComputer<storm::RationalNumber>;
+            template class DsMpiMdpUpperRewardBoundsComputer<storm::RationalNumber>;
+#endif
+        }
+    }
+}

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

@@ -0,0 +1,98 @@
+#pragma once
+
+#include <vector>
+
+namespace storm {
+    namespace storage {
+        template<typename ValueType>
+        class SparseMatrix;
+    }
+    
+    namespace modelchecker {
+        namespace helper {
+        
+            template<typename ValueType>
+            class DsMpiDtmcPriorityLess;
+            
+            template<typename ValueType>
+            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<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> const& oneStepTargetProbabilities);
+                
+                /*!
+                 * Computes upper bounds on the expected rewards.
+                 */
+                std::vector<ValueType> 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<ValueType> const& transitionMatrix;
+                std::vector<ValueType> const& originalRewards;
+                std::vector<ValueType> const& originalOneStepTargetProbabilities;
+                
+                // Derived from input data.
+                storm::storage::SparseMatrix<ValueType> backwardTransitions;
+                
+                // Data that the algorithm uses internally.
+                std::vector<ValueType> p;
+                std::vector<ValueType> w;
+                std::vector<ValueType> rewards;
+                std::vector<ValueType> targetProbabilities;
+                
+                friend class DsMpiDtmcPriorityLess<ValueType>;
+            };
+         
+            template<typename ValueType>
+            class DsMpiMdpPriorityLess;
+
+            template<typename ValueType>
+            class DsMpiMdpUpperRewardBoundsComputer : public DsMpiDtmcUpperRewardBoundsComputer<ValueType> {
+            public:
+                /*!
+                 * Creates an object that can compute upper bounds on the expected rewards for the provided DTMC.
+                 */
+                DsMpiMdpUpperRewardBoundsComputer(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& rewards, std::vector<ValueType> 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<uint64_t> choiceToState;
+                std::vector<uint64_t> policy;
+                
+                friend class DsMpiMdpPriorityLess<ValueType>;
+            };
+        }
+    }
+}

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<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), 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());
-                    auto one = storm::utility::one<ValueType>();
-                    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<ValueType>();
-                    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<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];
-                            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<ValueType>().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<PriorityLess> 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<ValueType> const& transitionMatrix;
-                std::vector<ValueType> const& originalRewards;
-                std::vector<ValueType> const& originalOneStepTargetProbabilities;
-                
-                // Derived from input data.
-                storm::storage::SparseMatrix<ValueType> backwardTransitions;
-                
-                // Data that the algorithm uses internally.
-                std::vector<ValueType> p;
-                std::vector<ValueType> w;
-                std::vector<ValueType> rewards;
-                std::vector<ValueType> targetProbabilities;
-            };
-            
             // 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);
+                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.");

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

@@ -2,9 +2,9 @@
 
 #include <boost/container/flat_map.hpp>
 
-
 #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<typename ValueType>
             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<bool>(upperResultBound);
                 }
+                
+                bool hasUpperResultBounds() const {
+                    return static_cast<bool>(upperResultBounds);
+                }
 
                 ValueType const& getUpperResultBound() const {
                     return upperResultBound.get();
                 }
+
+                std::vector<ValueType>& getUpperResultBounds() {
+                    return upperResultBounds.get();
+                }
+                
+                std::vector<ValueType> const& getUpperResultBounds() const {
+                    return upperResultBounds.get();
+                }
                 
                 std::vector<uint64_t>& getSchedulerHint() {
                     return schedulerHint.get();
@@ -150,11 +162,17 @@ namespace storm {
                     return eliminateEndComponents;
                 }
 
+                bool getComputeUpperBounds() {
+                    return computeUpperBounds;
+                }
+                
                 boost::optional<std::vector<uint64_t>> schedulerHint;
                 boost::optional<std::vector<ValueType>> valueHint;
                 boost::optional<ValueType> lowerResultBound;
                 boost::optional<ValueType> upperResultBound;
+                boost::optional<std::vector<ValueType>> upperResultBounds;
                 bool eliminateEndComponents;
+                bool computeUpperBounds;
             };
             
             template<typename ValueType>
@@ -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<ValueType>();
                 }
+                
+                // 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<ValueType>::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<ValueType> result(std::move(x));
 
@@ -485,23 +526,29 @@ namespace storm {
             };
             
             template<typename ValueType>
-            SparseMdpEndComponentInformation<ValueType> eliminateEndComponents(storm::storage::MaximalEndComponentDecomposition<ValueType> const& endComponentDecomposition, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const* sumColumns, storm::storage::BitVector const* selectedChoices, std::vector<ValueType> const* summand, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
+            SparseMdpEndComponentInformation<ValueType> eliminateEndComponents(storm::storage::MaximalEndComponentDecomposition<ValueType> const& endComponentDecomposition, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const* sumColumns, storm::storage::BitVector const* selectedChoices, std::vector<ValueType> const* summand, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>* columnSumVector, std::vector<ValueType>* summandResultVector) {
 
                 SparseMdpEndComponentInformation<ValueType> result(endComponentDecomposition, maybeStates);
 
                 // (1) Compute the number of maybe states not in ECs before any other maybe state.
                 std::vector<uint64_t> maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices();
-
-                // Create temporary vector storing possible transitions to ECs.
-                std::vector<std::pair<uint64_t, ValueType>> 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<ValueType> 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<std::pair<uint64_t, ValueType>> 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<ValueType>(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, b);
+                    return eliminateEndComponents<ValueType>(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<typename ValueType>
-            void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
+            void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, std::vector<ValueType>* 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<typename ValueType>
-            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
+            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, boost::optional<std::vector<ValueType>>& 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<ValueType>(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, b);
+                    return eliminateEndComponents<ValueType>(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<typename ValueType>
+            void computeUpperRewardBounds(SparseMdpHintType<ValueType>& hintInformation, storm::OptimizationDirection const& direction, storm::storage::SparseMatrix<ValueType> const& submatrix, std::vector<ValueType> const& choiceRewards, std::vector<ValueType> 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<ValueType> dsmpi(submatrix, choiceRewards, oneStepTargetProbabilities);
+                    hintInformation.upperResultBounds = dsmpi.computeUpperBounds();
+                } else {
+                    STORM_LOG_ASSERT(false, "Not yet implemented.");
+                }
+            }
+            
             template<typename ValueType>
             MDPSparseModelCheckingHelperReturnType<ValueType> SparseMdpPrctlHelper<ValueType>::computeReachabilityRewardsHelper(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint) {
                 
@@ -997,17 +1062,30 @@ namespace storm {
                         // Declare the components of the equation system we will solve.
                         storm::storage::SparseMatrix<ValueType> submatrix;
                         std::vector<ValueType> b;
+
+                        // If we need to compute upper bounds on the reward values, we need the one step probabilities
+                        // to a target state.
+                        boost::optional<std::vector<ValueType>> oneStepTargetProbabilities;
+                        if (hintInformation.getComputeUpperBounds()) {
+                            oneStepTargetProbabilities = std::vector<ValueType>();
+                        }
                         
                         // If the hint information tells us that we have to eliminate MECs, we do so now.
                         boost::optional<SparseMdpEndComponentInformation<ValueType>> 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.

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<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setTerminationCondition(std::unique_ptr<TerminationCondition<ValueType>> terminationCondition) {
+            this->terminationCondition = std::move(terminationCondition);
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::resetTerminationCondition() {
+            this->terminationCondition = nullptr;
+        }
+        
+        template<typename ValueType>
+        bool AbstractEquationSolver<ValueType>::hasCustomTerminationCondition() const {
+            return static_cast<bool>(this->terminationCondition);
+        }
+        
+        template<typename ValueType>
+        TerminationCondition<ValueType> const& AbstractEquationSolver<ValueType>::getTerminationCondition() const {
+            return *terminationCondition;
+        }
+        
+        template<typename ValueType>
+        bool AbstractEquationSolver<ValueType>::hasLowerBound(BoundType const& type) const {
+            if (type == BoundType::Any) {
+                return static_cast<bool>(lowerBound) || static_cast<bool>(lowerBounds);
+            } else if (type == BoundType::Global) {
+                return static_cast<bool>(lowerBound);
+            } else if (type == BoundType::Local) {
+                return static_cast<bool>(lowerBounds);
+            }
+            return false;
+        }
+        
+        template<typename ValueType>
+        bool AbstractEquationSolver<ValueType>::hasUpperBound(BoundType const& type) const {
+            if (type == BoundType::Any) {
+                return static_cast<bool>(upperBound) || static_cast<bool>(upperBounds);
+            } else if (type == BoundType::Global) {
+                return static_cast<bool>(upperBound);
+            } else if (type == BoundType::Local) {
+                return static_cast<bool>(upperBounds);
+            }
+            return false;
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setLowerBound(ValueType const& value) {
+            lowerBound = value;
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setUpperBound(ValueType const& value) {
+            upperBound = value;
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setBounds(ValueType const& lower, ValueType const& upper) {
+            setLowerBound(lower);
+            setUpperBound(upper);
+        }
+        
+        template<typename ValueType>
+        ValueType const& AbstractEquationSolver<ValueType>::getLowerBound() const {
+            return lowerBound.get();
+        }
+        
+        template<typename ValueType>
+        ValueType const& AbstractEquationSolver<ValueType>::getUpperBound() const {
+            return upperBound.get();
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType> const& AbstractEquationSolver<ValueType>::getLowerBounds() const {
+            return lowerBounds.get();
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType> const& AbstractEquationSolver<ValueType>::getUpperBounds() const {
+            return upperBounds.get();
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setLowerBounds(std::vector<ValueType> const& values) {
+            lowerBounds = values;
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType> const& values) {
+            upperBounds = values;
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType>&& values) {
+            upperBounds = std::move(values);
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setBounds(std::vector<ValueType> const& lower, std::vector<ValueType> const& upper) {
+            setLowerBounds(lower);
+            setUpperBounds(upper);
+        }
+        
+        template<typename ValueType>
+        void AbstractEquationSolver<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 AbstractEquationSolver<ValueType>::createUpperBoundsVector(std::unique_ptr<std::vector<ValueType>>& 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<std::vector<ValueType>>(this->getUpperBounds());
+                } else {
+                    upperBoundsVector = std::make_unique<std::vector<ValueType>>(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<double>;
+        
+#ifdef STORM_HAVE_CARL
+        template class AbstractEquationSolver<storm::RationalNumber>;
+        template class AbstractEquationSolver<storm::RationalFunction>;
+#endif
+        
+    }
+}

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 <memory>
 
+#include <boost/optional.hpp>
+
+#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<ValueType>> terminationCondition) {
-                this->terminationCondition = std::move(terminationCondition);
-            }
+            void setTerminationCondition(std::unique_ptr<TerminationCondition<ValueType>> 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<bool>(this->terminationCondition);
-            }
+            bool hasCustomTerminationCondition() const;
             
             /*!
              * Retrieves the custom termination condition (if any was set).
              * 
              * @return The custom termination condition.
              */
-            TerminationCondition<ValueType> const& getTerminationCondition() const {
-                return *terminationCondition;
-            }
+            TerminationCondition<ValueType> 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<ValueType> const& getLowerBounds() const;
+            
+            /*!
+             * Retrieves a vector containing the upper bounds (if there are any).
+             */
+            std::vector<ValueType> const& getUpperBounds() const;
+            
+            /*!
+             * Sets lower bounds for the solution that can potentially be used by the solver.
+             */
+            void setLowerBounds(std::vector<ValueType> const& values);
+            
+            /*!
+             * Sets upper bounds for the solution that can potentially be used by the solver.
+             */
+            void setUpperBounds(std::vector<ValueType> const& values);
+            
+            /*!
+             * Sets upper bounds for the solution that can potentially be used by the solver.
+             */
+            void setUpperBounds(std::vector<ValueType>&& values);
+            
+            /*!
+             * Sets bounds for the solution that can potentially be used by the solver.
+             */
+            void setBounds(std::vector<ValueType> const& lower, std::vector<ValueType> const& upper);
             
         protected:
+            void createUpperBoundsVector(std::unique_ptr<std::vector<ValueType>>& upperBoundsVector, uint64_t length) const;
+            void createLowerBoundsVector(std::vector<ValueType>& lowerBoundsVector) const;
+
             // A termination condition to be used (can be unset).
             std::unique_ptr<TerminationCondition<ValueType>> terminationCondition;
+            
+            // A lower bound if one was set.
+            boost::optional<ValueType> lowerBound;
+            
+            // An upper bound if one was set.
+            boost::optional<ValueType> upperBound;
+            
+            // Lower bounds if they were set.
+            boost::optional<std::vector<ValueType>> lowerBounds;
+            
+            // Lower bounds if they were set.
+            boost::optional<std::vector<ValueType>> upperBounds;
         };
         
     }

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -380,37 +380,12 @@ namespace storm {
                 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(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<ValueType> submatrix = this->A->selectRowsFromRowGroups(this->getInitialScheduler(), convertToEquationSystem);
-                if (convertToEquationSystem) {
-                    submatrix.convertToEquationSystem();
-                }
-                storm::utility::vector::selectVectorValues<ValueType>(*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<ValueType>* lowerX = &x;
+            this->createLowerBoundsVector(*lowerX);
+            this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount());
             std::vector<ValueType>* upperX = this->auxiliaryRowGroupVector.get();
             std::vector<ValueType>* tmp = nullptr;
             if (!useGaussSeidelMultiplication) {

src/storm/solver/LinearEquationSolver.cpp

@@ -132,122 +132,6 @@ namespace storm {
             cachedRowVector.reset();
         }
         
-        template<typename ValueType>
-        bool LinearEquationSolver<ValueType>::hasLowerBound(BoundType const& type) const {
-            if (type == BoundType::Any) {
-                return static_cast<bool>(lowerBound) || static_cast<bool>(lowerBounds);
-            } else if (type == BoundType::Global) {
-                return static_cast<bool>(lowerBound);
-            } else if (type == BoundType::Local) {
-                return static_cast<bool>(lowerBounds);
-            }
-            return false;
-        }
-        
-        template<typename ValueType>
-        bool LinearEquationSolver<ValueType>::hasUpperBound(BoundType const& type) const {
-            if (type == BoundType::Any) {
-                return static_cast<bool>(upperBound) || static_cast<bool>(upperBounds);
-            } else if (type == BoundType::Global) {
-                return static_cast<bool>(upperBound);
-            } else if (type == BoundType::Local) {
-                return static_cast<bool>(upperBounds);
-            }
-            return false;
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setLowerBound(ValueType const& value) {
-            lowerBound = value;
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setUpperBound(ValueType const& value) {
-            upperBound = value;
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setBounds(ValueType const& lower, ValueType const& upper) {
-            setLowerBound(lower);
-            setUpperBound(upper);
-        }
-        
-        template<typename ValueType>
-        ValueType const& LinearEquationSolver<ValueType>::getLowerBound() const {
-            return lowerBound.get();
-        }
-        
-        template<typename ValueType>
-        ValueType const& LinearEquationSolver<ValueType>::getUpperBound() const {
-            return upperBound.get();
-        }
-        
-        template<typename ValueType>
-        std::vector<ValueType> const& LinearEquationSolver<ValueType>::getLowerBounds() const {
-            return lowerBounds.get();
-        }
-        
-        template<typename ValueType>
-        std::vector<ValueType> const& LinearEquationSolver<ValueType>::getUpperBounds() const {
-            return upperBounds.get();
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setLowerBounds(std::vector<ValueType> const& values) {
-            lowerBounds = values;
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType> const& values) {
-            upperBounds = values;
-        }
-
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType>&& values) {
-            upperBounds = std::move(values);
-        }
-
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::setBounds(std::vector<ValueType> const& lower, std::vector<ValueType> const& upper) {
-            setLowerBounds(lower);
-            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) {
-                if (this->hasUpperBound(BoundType::Local)) {
-                    upperBoundsVector = std::make_unique<std::vector<ValueType>>(this->getUpperBounds());
-                } else {
-                    upperBoundsVector = std::make_unique<std::vector<ValueType>>(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<typename ValueType>
         std::unique_ptr<LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
             std::unique_ptr<LinearEquationSolver<ValueType>> solver = this->create();

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<ValueType> const& getLowerBounds() const;
-            
-            /*!
-             * Retrieves a vector containing the upper bounds (if there are any).
-             */
-            std::vector<ValueType> const& getUpperBounds() const;
-            
-            /*!
-             * Sets lower bounds for the solution that can potentially be used by the solver.
-             */
-            void setLowerBounds(std::vector<ValueType> const& values);
-            
-            /*!
-             * Sets upper bounds for the solution that can potentially be used by the solver.
-             */
-            void setUpperBounds(std::vector<ValueType> const& values);
-
-            /*!
-             * Sets upper bounds for the solution that can potentially be used by the solver.
-             */
-            void setUpperBounds(std::vector<ValueType>&& values);
-
-            /*!
-             * Sets bounds for the solution that can potentially be used by the solver.
-             */
-            void setBounds(std::vector<ValueType> const& lower, std::vector<ValueType> const& upper);
-            
         protected:
             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;
             
-            // A lower bound if one was set.
-            boost::optional<ValueType> lowerBound;
-
-            // An upper bound if one was set.
-            boost::optional<ValueType> upperBound;
-            
-            // Lower bounds if they were set.
-            boost::optional<std::vector<ValueType>> lowerBounds;
-
-            // Lower bounds if they were set.
-            boost::optional<std::vector<ValueType>> upperBounds;
-
         private:
             /*!
              * Retrieves the row count of the matrix associated with this solver.

src/storm/solver/MinMaxLinearEquationSolver.cpp

@@ -110,43 +110,7 @@ namespace storm {
         void MinMaxLinearEquationSolver<ValueType>::clearCache() const {
             // Intentionally left empty.
         }
-        
-        template<typename ValueType>
-        void MinMaxLinearEquationSolver<ValueType>::setLowerBound(ValueType const& value) {
-            lowerBound = value;
-        }
-        
-        template<typename ValueType>
-        void MinMaxLinearEquationSolver<ValueType>::setUpperBound(ValueType const& value) {
-            upperBound = value;
-        }
-        
-        template<typename ValueType>
-        void MinMaxLinearEquationSolver<ValueType>::setBounds(ValueType const& lower, ValueType const& upper) {
-            setLowerBound(lower);
-            setUpperBound(upper);
-        }
-        
-        template<typename ValueType>
-        bool MinMaxLinearEquationSolver<ValueType>::hasUpperBound() const {
-            return static_cast<bool>(upperBound);
-        }
-        
-        template<typename ValueType>
-        bool MinMaxLinearEquationSolver<ValueType>::hasLowerBound() const {
-            return static_cast<bool>(lowerBound);
-        }
-
-        template<typename ValueType>
-        ValueType const& MinMaxLinearEquationSolver<ValueType>::getUpperBound() const {
-            return upperBound.get();
-        }
-        
-        template<typename ValueType>
-        ValueType const& MinMaxLinearEquationSolver<ValueType>::getLowerBound() const {
-            return lowerBound.get();
-        }
-        
+                
         template<typename ValueType>
         void MinMaxLinearEquationSolver<ValueType>::setInitialScheduler(std::vector<uint_fast64_t>&& choices) {
             initialScheduler = std::move(choices);

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<ValueType>& x, std::vector<ValueType> 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<std::vector<uint_fast64_t>> schedulerChoices;
             
-            // A lower bound if one was set.
-            boost::optional<ValueType> lowerBound;
-            
-            // An upper bound if one was set.
-            boost::optional<ValueType> upperBound;
-            
             // A scheduler that can be used by solvers that require a valid initial scheduler.
             boost::optional<std::vector<uint_fast64_t>> initialScheduler;
             

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -438,7 +438,7 @@ namespace storm {
             
             std::vector<ValueType>* lowerX = &x;
             this->createLowerBoundsVector(*lowerX);
-            this->createUpperBoundsVector(this->cachedRowVector);
+            this->createUpperBoundsVector(this->cachedRowVector, this->getMatrixRowCount());
             std::vector<ValueType>* upperX = this->cachedRowVector.get();
             
             bool useGaussSeidelMultiplication = this->getSettings().getPowerMethodMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;

src/storm/storage/ConsecutiveUint64DynamicPriorityQueue.h

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <vector>
+#include <numeric>
 
 #include "storm/utility/macros.h"