respecting lower/upper bounds from preprocessing in quick sound power method

src/storm/solver/AbstractEquationSolver.cpp

@@ -9,6 +9,7 @@
 #include "storm/utility/constants.h"
 #include "storm/utility/macros.h"
 #include "storm/exceptions/UnmetRequirementException.h"
+#include "storm/exceptions/InvalidOperationException.h"
 
 namespace storm {
     namespace solver {
@@ -119,10 +120,33 @@ namespace storm {
             return lowerBound.get();
         }
         
+        template<typename ValueType>
+        ValueType AbstractEquationSolver<ValueType>::getLowerBound(bool convertLocalBounds) const {
+            if (lowerBound) {
+                return lowerBound.get();
+            } else if (convertLocalBounds) {
+                return *std::min_element(lowerBounds->begin(), lowerBounds->end());
+            }
+            STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "No lower bound available but some was requested.");
+            return ValueType();
+        }
+        
         template<typename ValueType>
         ValueType const& AbstractEquationSolver<ValueType>::getUpperBound() const {
             return upperBound.get();
         }
+
+        template<typename ValueType>
+        ValueType AbstractEquationSolver<ValueType>::getUpperBound(bool convertLocalBounds) const {
+            if (upperBound) {
+                return upperBound.get();
+            } else if (convertLocalBounds) {
+                return *std::max_element(upperBounds->begin(), upperBounds->end());
+            }
+            STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "No upper bound available but some was requested.");
+            return ValueType();
+        }
+
         
         template<typename ValueType>
         std::vector<ValueType> const& AbstractEquationSolver<ValueType>::getLowerBounds() const {

src/storm/solver/AbstractEquationSolver.h

@@ -97,11 +97,25 @@ namespace storm {
              */
             ValueType const& getLowerBound() const;
             
+            /*!
+             * Retrieves the lower bound (if there is any).
+             * If the given flag is true and if there are only local bounds,
+             * the minimum of the local bounds is returned.
+             */
+            ValueType getLowerBound(bool convertLocalBounds) const;
+            
             /*!
              * Retrieves the upper bound (if there is any).
              */
             ValueType const& getUpperBound() const;
             
+           /*!
+             * Retrieves the upper bound (if there is any).
+             * If the given flag is true and if there are only local bounds,
+             * the maximum of the local bounds is returned.
+             */
+            ValueType getUpperBound(bool convertLocalBounds) const;
+            
             /*!
              * Retrieves a vector containing the lower bounds (if there are any).
              */

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -984,20 +984,13 @@ namespace storm {
             
             
             // Prepare initial bounds for the solution (if given)
-            if (this->hasLowerBound(AbstractEquationSolver<ValueType>::BoundType::Global)) {
-                helper.setLowerBound(this->getLowerBound());
-            } else if (this->hasLowerBound(AbstractEquationSolver<ValueType>::BoundType::Local)) {
-                helper.setLowerBound(*std::min_element(this->getLowerBounds().begin(), this->getLowerBounds().end()));
+            if (this->hasLowerBound()) {
+                helper.setLowerBound(this->getLowerBound(true));
             }
-            if (this->hasUpperBound(AbstractEquationSolver<ValueType>::BoundType::Global)) {
-                helper.setUpperBound(this->getUpperBound());
-            } else if (this->hasUpperBound(AbstractEquationSolver<ValueType>::BoundType::Local)) {
-                helper.setUpperBound(*std::max_element(this->getUpperBounds().begin(), this->getUpperBounds().end()));
+            if (this->hasUpperBound()) {
+                helper.setUpperBound(this->getUpperBound(true));
             }
             
-            //STORM_LOG_INFO_COND(!hasCurrentLowerBound, "Initial lower bound on the result is " << currentLowerBound);
-            //STORM_LOG_INFO_COND(!hasCurrentUpperBound, "Initial upper bound on the result is " << currentUpperBound);
-            
             storm::storage::BitVector const* relevantValuesPtr = nullptr;
             if (this->hasRelevantValues()) {
                 relevantValuesPtr = &this->getRelevantValues();

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -604,8 +604,19 @@ namespace storm {
             uint64_t iterations = 0;
             bool converged = false;
             bool terminate = false;
-            ValueType minValueBound, maxValueBound;
             uint64_t minIndex(0), maxIndex(0);
+            ValueType minValueBound, maxValueBound;
+            bool hasMinValueBound, hasMaxValueBound;
+            // Prepare initial bounds for the solution (if given)
+            if (this->hasLowerBound()) {
+                minValueBound = this->getLowerBound(true);
+                hasMinValueBound = true;
+            }
+            if (this->hasUpperBound()) {
+                maxValueBound = this->getUpperBound(true);
+                hasMaxValueBound = true;
+            }
+            
             bool convergencePhase1 = true;
             uint64_t firstIndexViolatingConvergence = 0;
             this->startMeasureProgress();
@@ -646,11 +657,17 @@ namespace storm {
                 if (!convergencePhase1) {
                     // Phase 2: the difference between lower and upper bound has to be < precision at every (relevant) value
                     // First check with (possibly too tight) bounds from a previous iteration. Only compute the actual bounds if this first check passes.
-                    minValueBound = stepBoundedX->at(minIndex) / (storm::utility::one<ValueType>() - stepBoundedStayProbs->at(minIndex));
-                    maxValueBound = stepBoundedX->at(maxIndex) / (storm::utility::one<ValueType>() - stepBoundedStayProbs->at(maxIndex));
+                    ValueType minValueBoundCandidate = stepBoundedX->at(minIndex) / (storm::utility::one<ValueType>() - stepBoundedStayProbs->at(minIndex));
+                    ValueType maxValueBoundCandidate = stepBoundedX->at(maxIndex) / (storm::utility::one<ValueType>() - stepBoundedStayProbs->at(maxIndex));
+                    if (hasMinValueBound && minValueBound > minValueBoundCandidate) {
+                        minValueBoundCandidate = minValueBound;
+                    }
+                    if (hasMaxValueBound && maxValueBound < maxValueBoundCandidate) {
+                        maxValueBoundCandidate = maxValueBound;
+                    }
                     ValueType const& stayProb = stepBoundedStayProbs->at(firstIndexViolatingConvergence);
                     // The error made in this iteration
-                    ValueType absoluteError = stayProb * (maxValueBound - minValueBound);
+                    ValueType absoluteError = stayProb * (maxValueBoundCandidate - minValueBoundCandidate);
                     // The maximal allowed error (possibly respecting relative precision)
                     // Note: We implement the relative convergence criterion in a way that avoids division by zero in the case where stepBoundedX[i] is zero.
                     ValueType maxAllowedError = relative ? (precision * stepBoundedX->at(firstIndexViolatingConvergence)) : precision;
@@ -661,14 +678,20 @@ namespace storm {
                         auto probIt = stepBoundedStayProbs->begin();
                         for (uint64_t index = 0; valIt != valIte; ++valIt, ++probIt, ++index) {
                             ValueType currentBound = *valIt / (storm::utility::one<ValueType>() - *probIt);
-                            if (currentBound < minValueBound) {
+                            if (currentBound < minValueBoundCandidate) {
                                 minIndex = index;
-                                minValueBound = std::move(currentBound);
-                            } else if (currentBound > maxValueBound) {
+                                minValueBoundCandidate = std::move(currentBound);
+                            } else if (currentBound > maxValueBoundCandidate) {
                                 maxIndex = index;
-                                maxValueBound = std::move(currentBound);
+                                maxValueBoundCandidate = std::move(currentBound);
                             }
                         }
+                        if (!hasMinValueBound || minValueBoundCandidate > minValueBound) {
+                            minValueBound = minValueBoundCandidate;
+                        }
+                        if (!hasMaxValueBound || maxValueBoundCandidate < maxValueBound) {
+                            maxValueBound = maxValueBoundCandidate;
+                        }
                         absoluteError = stayProb * (maxValueBound - minValueBound);
                         if (absoluteError <= maxAllowedError) {
                             // The current index satisfies the desired bound. We now move to the next index that violates it