Merge branch 'master' into monitoring

CHANGELOG.md

@@ -8,7 +8,7 @@ The releases of major and minor versions contain an overview of changes since th
 Version 1.6.x
 -------------
 
-## Version 1.6.1 (??)
+## Version 1.6.2 (2020/09)
 - Prism program simplification improved.
 - Revamped implementation of long-run-average algorithms, including scheduler export for LRA properties on Markov automata.
 - Support for step-bounded properties of the form ... [F[x,y] ... ] for DTMCs and MDPs (sparse engine).

CMakeLists.txt

@@ -61,6 +61,9 @@ export_option(STORM_USE_CLN_RF)
 option(BUILD_SHARED_LIBS "Build the Storm library dynamically" OFF)
 option(STORM_DEBUG_CUDD "Build CUDD in debug mode." OFF)
 MARK_AS_ADVANCED(STORM_DEBUG_CUDD)
+option(STORM_EXCLUDE_TESTS_FROM_ALL "If set, tests will not be compiled by default" OFF )
+export_option(STORM_EXCLUDE_TESTS_FROM_ALL)
+MARK_AS_ADVANCED(STORM_EXCLUDE_TESTS_FROM_ALL)
 set(BOOST_ROOT "" CACHE STRING "A hint to the root directory of Boost (optional).")
 set(GUROBI_ROOT "" CACHE STRING "A hint to the root directory of Gurobi (optional).")
 set(Z3_ROOT "" CACHE STRING "A hint to the root directory of Z3 (optional).")

doc/checklist_new_release.md

@@ -34,11 +34,18 @@ Note that in most cases a simultaneous release of [carl](https://github.com/smtr
 6. [Add new release](https://github.com/moves-rwth/storm/releases/new) in GitHub.
 
 7. Update `stable` branch:
+
    ```console
    git checkout stable
-   git merge master
+   git rebase master
    git push origin stable
    ```
+   Note: Rebasing might fail if `stable` is ahead of `master` (e.g. because of merge commits). In this case we can do:
+    ```console
+   git checkout stable
+   git reset --hard master
+   git push --force origin stable
+   ```
 
 8. Update [Homebrew formula](https://github.com/moves-rwth/homebrew-storm).
 

resources/3rdparty/CMakeLists.txt

@@ -260,14 +260,26 @@ if (NOT STORM_FORCE_SHIPPED_CARL)
 		find_package(carl QUIET)
     endif()
 endif()
+
+set(STORM_SHIPPED_CARL OFF)
 if(carl_FOUND AND NOT STORM_FORCE_SHIPPED_CARL)
     get_target_property(carlLOCATION lib_carl LOCATION)
     if("${carlLOCATION}" STREQUAL "carlLOCATION-NOTFOUND")
-        message(FATAL_ERROR "Library location for carl is not found, did you build carl?")
+        if(EXISTS ${STORM_3RDPARTY_BINARY_DIR}/carl)
+            message(WARNING "Storm - Library for carl location is not found but the directory ${STORM_3RDPARTY_BINARY_DIR}/carl exists. Will (re-)try to build a shipped version of carl.")
+            set(STORM_SHIPPED_CARL ON)
+        else()
+            message(FATAL_ERROR "Library location for carl is not found, did you build carl?")
+        endif()
     elseif(EXISTS ${carlLOCATION})
         #empty on purpose
     else()
-        message(FATAL_ERROR "File ${carlLOCATION} does not exist, did you build carl?")
+        if(EXISTS ${STORM_3RDPARTY_BINARY_DIR}/carl)
+            message(WARNING "Storm - File ${carlLOCATION} does not exist but the directory ${STORM_3RDPARTY_BINARY_DIR}/carl exists. Will (re-)try to build a shipped version of carl.")
+            set(STORM_SHIPPED_CARL ON)
+        else()
+            message(FATAL_ERROR "File ${carlLOCATION} does not exist, did you build carl?")
+        endif()
     endif()
     if("${carl_VERSION_MAJOR}" STREQUAL "${CARL_C14VERSION}")
         message(STATUS "Storm - Found carl using master14 branch.")
@@ -276,14 +288,16 @@ if(carl_FOUND AND NOT STORM_FORCE_SHIPPED_CARL)
         message(FATAL_ERROR "Carl outdated, require ${CARL_MINVERSION}, have ${carl_VERSION}")
     endif()
 
-    set(STORM_SHIPPED_CARL OFF)
     set(STORM_HAVE_CARL ON)
     message(STATUS "Storm - Use system version of carl.")
     message(STATUS "Storm - Linking with preinstalled carl ${carl_VERSION} (include: ${carl_INCLUDE_DIR}, library ${carl_LIBRARIES}, CARL_USE_CLN_NUMBERS: ${CARL_USE_CLN_NUMBERS}, CARL_USE_GINAC: ${CARL_USE_GINAC}).")
     set(STORM_HAVE_CLN ${CARL_USE_CLN_NUMBERS})
     set(STORM_HAVE_GINAC ${CARL_USE_GINAC})
 else()
-	set(STORM_SHIPPED_CARL ON)
+    set(STORM_SHIPPED_CARL ON)
+endif()
+
+if (STORM_SHIPPED_CARL)
 	# The first external project will be built at *configure stage*
 	message(STATUS "Carl - Start of config process")
 	file(MAKE_DIRECTORY ${STORM_3RDPARTY_BINARY_DIR}/carl_download)
@@ -333,7 +347,11 @@ else()
     message(STATUS "Storm - Linking with shipped carl ${carl_VERSION} (include: ${carl_INCLUDE_DIR}, library ${carl_LIBRARIES}, CARL_USE_CLN_NUMBERS: ${CARL_USE_CLN_NUMBERS}, CARL_USE_GINAC: ${CARL_USE_GINAC}).")
 
     # install the carl dynamic library if we built it
-    install(FILES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl.${carl_VERSION}${DYNAMIC_EXT} DESTINATION lib)
+    if(MACOSX)
+        install(FILES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl.${carl_VERSION}${DYNAMIC_EXT} DESTINATION lib)
+    else()
+        install(FILES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl${DYNAMIC_EXT}.${carl_VERSION} DESTINATION lib)
+    endif()
 endif()
 
 if(STORM_USE_CLN_RF AND NOT STORM_HAVE_CLN)

src/CMakeLists.txt

@@ -23,7 +23,10 @@ add_subdirectory(storm-pomdp-cli)
 add_subdirectory(storm-conv)
 add_subdirectory(storm-conv-cli)
 
-
-add_subdirectory(test)
+if (STORM_EXCLUDE_TESTS_FROM_ALL)
+    add_subdirectory(test EXCLUDE_FROM_ALL)
+else()
+    add_subdirectory(test)
+endif()
 
 set(STORM_TARGETS ${STORM_TARGETS} PARENT_SCOPE)

src/storm-pars/modelchecker/instantiation/SparseDtmcInstantiationModelChecker.cpp

@@ -42,8 +42,29 @@ namespace storm {
                 this->currentCheckTask->setHint(std::make_shared<ExplicitModelCheckerHint<ConstantType>>());
             }
             ExplicitModelCheckerHint<ConstantType>& hint = this->currentCheckTask->getHint().template asExplicitModelCheckerHint<ConstantType>();
-            std::unique_ptr<CheckResult> result;
             
+            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
+                // Perform purely qualitative analysis once
+                std::vector<ConstantType> qualitativeResult;
+                if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
+                    auto newCheckTask = *this->currentCheckTask;
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    qualitativeResult = modelChecker.check(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                } else {
+                    auto newCheckTask = this->currentCheckTask->substituteFormula(this->currentCheckTask->getFormula().asOperatorFormula().getSubformula());
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    qualitativeResult = modelChecker.computeProbabilities(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                }
+                storm::storage::BitVector maybeStates = storm::utility::vector::filter<ConstantType>(qualitativeResult,
+                                [] (ConstantType const& value) -> bool { return !(storm::utility::isZero<ConstantType>(value) || storm::utility::isOne<ConstantType>(value)); });
+                hint.setMaybeStates(std::move(maybeStates));
+                hint.setResultHint(std::move(qualitativeResult));
+                hint.setComputeOnlyMaybeStates(true);
+            }
+            
+            std::unique_ptr<CheckResult> result;
             // Check the formula and store the result as a hint for the next call.
             // For qualitative properties, we still want a quantitative result hint. Hence we perform the check on the subformula
             if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
@@ -56,15 +77,6 @@ namespace storm {
                 hint.setResultHint(std::move(quantitativeResult->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector()));
             }
             
-            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
-                // Extract the maybe states from the current result.
-                assert(hint.hasResultHint());
-                storm::storage::BitVector maybeStates = ~storm::utility::vector::filter<ConstantType>(hint.getResultHint(),
-                                [] (ConstantType const& value) -> bool { return storm::utility::isZero<ConstantType>(value) || storm::utility::isOne<ConstantType>(value); });
-                hint.setMaybeStates(std::move(maybeStates));
-                hint.setComputeOnlyMaybeStates(true);
-            }
-            
             return result;
         }
         
@@ -75,6 +87,28 @@ namespace storm {
                 this->currentCheckTask->setHint(std::make_shared<ExplicitModelCheckerHint<ConstantType>>());
             }
             ExplicitModelCheckerHint<ConstantType>& hint = this->currentCheckTask->getHint().template asExplicitModelCheckerHint<ConstantType>();
+            
+            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
+                // Perform purely qualitative analysis once
+                std::vector<ConstantType> qualitativeResult;
+                if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
+                    auto newCheckTask = *this->currentCheckTask;
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    qualitativeResult = modelChecker.check(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                } else {
+                    auto newCheckTask = this->currentCheckTask->substituteFormula(this->currentCheckTask->getFormula().asOperatorFormula().getSubformula());
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    qualitativeResult = modelChecker.computeRewards(env, this->currentCheckTask->getFormula().asRewardOperatorFormula().getMeasureType(), newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                }
+                storm::storage::BitVector maybeStates = storm::utility::vector::filter<ConstantType>(qualitativeResult,
+                                [] (ConstantType const& value) -> bool { return !(storm::utility::isZero<ConstantType>(value) || storm::utility::isInfinity<ConstantType>(value)); });
+                hint.setMaybeStates(std::move(maybeStates));
+                hint.setResultHint(std::move(qualitativeResult));
+                hint.setComputeOnlyMaybeStates(true);
+            }
+            
             std::unique_ptr<CheckResult> result;
             
             // Check the formula and store the result as a hint for the next call.
@@ -89,18 +123,6 @@ namespace storm {
                 this->currentCheckTask->getHint().template asExplicitModelCheckerHint<ConstantType>().setResultHint(std::move(quantitativeResult->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector()));
             }
             
-            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
-                // Extract the maybe states from the current result.
-                assert(hint.hasResultHint());
-                storm::storage::BitVector maybeStates = ~storm::utility::vector::filterInfinity(hint.getResultHint());
-                // We need to exclude the target states from the maybe states.
-                // Note that we can not consider the states with reward zero since a valuation might set a reward to zero
-                std::unique_ptr<CheckResult> subFormulaResult = modelChecker.check(env, this->currentCheckTask->getFormula().asOperatorFormula().getSubformula().asEventuallyFormula().getSubformula());
-                maybeStates = maybeStates & ~(subFormulaResult->asExplicitQualitativeCheckResult().getTruthValuesVector());
-                hint.setMaybeStates(std::move(maybeStates));
-                hint.setComputeOnlyMaybeStates(true);
-            }
-            
             return result;
         }
         

src/storm-pars/modelchecker/instantiation/SparseMdpInstantiationModelChecker.cpp

@@ -47,8 +47,31 @@ namespace storm {
                 this->currentCheckTask->setHint(std::make_shared<ExplicitModelCheckerHint<ConstantType>>());
             }
             ExplicitModelCheckerHint<ConstantType>& hint = this->currentCheckTask->getHint().template asExplicitModelCheckerHint<ConstantType>();
-            std::unique_ptr<CheckResult> result;
             
+            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
+                // Perform purely qualitative analysis once
+                std::vector<ConstantType> qualitativeResult;
+                if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
+                    auto newCheckTask = *this->currentCheckTask;
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    newCheckTask.setProduceSchedulers(false);
+                    qualitativeResult = modelChecker.check(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                } else {
+                    auto newCheckTask = this->currentCheckTask->substituteFormula(this->currentCheckTask->getFormula().asOperatorFormula().getSubformula());
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    newCheckTask.setProduceSchedulers(false);
+                    qualitativeResult = modelChecker.computeProbabilities(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                }
+                storm::storage::BitVector maybeStates = storm::utility::vector::filter<ConstantType>(qualitativeResult,
+                                [] (ConstantType const& value) -> bool { return !(storm::utility::isZero<ConstantType>(value) || storm::utility::isOne<ConstantType>(value)); });
+                hint.setMaybeStates(std::move(maybeStates));
+                hint.setResultHint(std::move(qualitativeResult));
+                hint.setComputeOnlyMaybeStates(true);
+            }
+            
+            std::unique_ptr<CheckResult> result;
             // Check the formula and store the result as a hint for the next call.
             // For qualitative properties, we still want a quantitative result hint. Hence we perform the check on the subformula
             if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
@@ -65,21 +88,6 @@ namespace storm {
                 hint.setSchedulerHint(std::move(dynamic_cast<storm::storage::Scheduler<ConstantType>&>(scheduler)));
             }
             
-            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
-                // Extract the maybe states from the current result.
-                assert(hint.hasResultHint());
-                storm::storage::BitVector maybeStates = ~storm::utility::vector::filter<ConstantType>(hint.getResultHint(),
-                                [] (ConstantType const& value) -> bool { return storm::utility::isZero<ConstantType>(value) || storm::utility::isOne<ConstantType>(value); });
-                hint.setMaybeStates(std::move(maybeStates));
-                hint.setComputeOnlyMaybeStates(true);
-                
-                // Check if there can be end components within the maybestates
-                if (storm::solver::minimize(this->currentCheckTask->getOptimizationDirection()) ||
-                storm::utility::graph::performProb1A(instantiatedModel.getTransitionMatrix(), instantiatedModel.getTransitionMatrix().getRowGroupIndices(), instantiatedModel.getBackwardTransitions(), hint.getMaybeStates(), ~hint.getMaybeStates()).full()) {
-                    hint.setNoEndComponentsInMaybeStates(true);
-                }
-            }
-            
             return result;
         }
         
@@ -92,8 +100,31 @@ namespace storm {
                 this->currentCheckTask->setHint(std::make_shared<ExplicitModelCheckerHint<ConstantType>>());
             }
             ExplicitModelCheckerHint<ConstantType>& hint = this->currentCheckTask->getHint().template asExplicitModelCheckerHint<ConstantType>();
-            std::unique_ptr<CheckResult> result;
             
+            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
+                // Perform purely qualitative analysis once
+                std::vector<ConstantType> qualitativeResult;
+                if (this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
+                    auto newCheckTask = *this->currentCheckTask;
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    newCheckTask.setProduceSchedulers(false);
+                    qualitativeResult = modelChecker.check(env, newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                } else {
+                    auto newCheckTask = this->currentCheckTask->substituteFormula(this->currentCheckTask->getFormula().asOperatorFormula().getSubformula());
+                    newCheckTask.setQualitative(true);
+                    newCheckTask.setOnlyInitialStatesRelevant(false);
+                    newCheckTask.setProduceSchedulers(false);
+                    qualitativeResult = modelChecker.computeRewards(env, this->currentCheckTask->getFormula().asRewardOperatorFormula().getMeasureType(), newCheckTask)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
+                }
+                storm::storage::BitVector maybeStates = storm::utility::vector::filter<ConstantType>(qualitativeResult,
+                                [] (ConstantType const& value) -> bool { return !(storm::utility::isZero<ConstantType>(value) || storm::utility::isInfinity<ConstantType>(value)); });
+                hint.setMaybeStates(std::move(maybeStates));
+                hint.setResultHint(std::move(qualitativeResult));
+                hint.setComputeOnlyMaybeStates(true);
+            }
+            
+            std::unique_ptr<CheckResult> result;
             // Check the formula and store the result as a hint for the next call.
             // For qualitative properties, we still want a quantitative result hint. Hence we perform the check on the subformula
             if(this->currentCheckTask->getFormula().asOperatorFormula().hasQuantitativeResult()) {
@@ -110,23 +141,6 @@ namespace storm {
                 hint.setSchedulerHint(std::move(dynamic_cast<storm::storage::Scheduler<ConstantType>&>(scheduler)));
             }
             
-            if (this->getInstantiationsAreGraphPreserving() && !hint.hasMaybeStates()) {
-                // Extract the maybe states from the current result.
-                assert(hint.hasResultHint());
-                storm::storage::BitVector maybeStates = ~storm::utility::vector::filterInfinity(hint.getResultHint());
-                // We need to exclude the target states from the maybe states.
-                // Note that we can not consider the states with reward zero since a valuation might set a reward to zero
-                std::unique_ptr<CheckResult> subFormulaResult = modelChecker.check(env, this->currentCheckTask->getFormula().asOperatorFormula().getSubformula().asEventuallyFormula().getSubformula());
-                maybeStates = maybeStates & ~(subFormulaResult->asExplicitQualitativeCheckResult().getTruthValuesVector());
-                hint.setMaybeStates(std::move(maybeStates));
-                hint.setComputeOnlyMaybeStates(true);
-                
-                // Check if there can be end components within the maybestates
-                if (storm::solver::maximize(this->currentCheckTask->getOptimizationDirection()) ||
-                storm::utility::graph::performProb1A(instantiatedModel.getTransitionMatrix(), instantiatedModel.getTransitionMatrix().getRowGroupIndices(), instantiatedModel.getBackwardTransitions(), hint.getMaybeStates(), ~hint.getMaybeStates()).full()) {
-                    hint.setNoEndComponentsInMaybeStates(true);
-                }
-            }
             return result;
         }
         

src/storm-pars/storage/ParameterRegion.cpp

@@ -1,7 +1,10 @@
 #include "storm-pars/storage/ParameterRegion.h"
 
+#include <limits>
+
 #include "storm/utility/macros.h"
 #include "storm/exceptions/InvalidArgumentException.h"
+#include "storm/exceptions/OutOfRangeException.h"
 #include "storm/utility/constants.h"
 
 namespace storm {
@@ -88,6 +91,7 @@ namespace storm {
         template<typename ParametricType>
         std::vector<typename ParameterRegion<ParametricType>::Valuation> ParameterRegion<ParametricType>::getVerticesOfRegion(std::set<VariableType> const& consideredVariables) const {
             std::size_t const numOfVariables = consideredVariables.size();
+            STORM_LOG_THROW(numOfVariables <= std::numeric_limits<std::size_t>::digits, storm::exceptions::OutOfRangeException, "Number of variables " << numOfVariables << " is too high.");
             std::size_t const numOfVertices = std::pow(2, numOfVariables);
             std::vector<Valuation> resultingVector(numOfVertices);
             

src/storm/environment/solver/OviSolverEnvironment.cpp

@@ -11,7 +11,6 @@ namespace storm {
         auto const& oviSettings = storm::settings::getModule<storm::settings::modules::OviSolverSettings>();
         precisionUpdateFactor = storm::utility::convertNumber<storm::RationalNumber>(oviSettings.getPrecisionUpdateFactor());
         maxVerificationIterationFactor = storm::utility::convertNumber<storm::RationalNumber>(oviSettings.getMaxVerificationIterationFactor());
-        relevantValuesForPrecisionUpdate = oviSettings.useRelevantValuesForPrecisionUpdate();
         upperBoundGuessingFactor = storm::utility::convertNumber<storm::RationalNumber>(oviSettings.getUpperBoundGuessingFactor());
         upperBoundOnlyIterations = oviSettings.getUpperBoundOnlyIterations();
         noTerminationGuaranteeMinimumMethod = oviSettings.useNoTerminationGuaranteeMinimumMethod();
@@ -29,10 +28,6 @@ namespace storm {
         return maxVerificationIterationFactor;
     }
     
-    bool OviSolverEnvironment::useRelevantValuesForPrecisionUpdate() const {
-        return relevantValuesForPrecisionUpdate;
-    }
-
     storm::RationalNumber OviSolverEnvironment::getUpperBoundGuessingFactor() const {
         return upperBoundGuessingFactor;
     }

src/storm/environment/solver/OviSolverEnvironment.h

@@ -14,7 +14,6 @@ namespace storm {
         
         storm::RationalNumber getPrecisionUpdateFactor() const;
         storm::RationalNumber getMaxVerificationIterationFactor() const;
-        bool useRelevantValuesForPrecisionUpdate() const;
         storm::RationalNumber getUpperBoundGuessingFactor() const;
         uint64_t getUpperBoundOnlyIterations() const;
         bool useNoTerminationGuaranteeMinimumMethod() const;
@@ -22,7 +21,6 @@ namespace storm {
     private:
         storm::RationalNumber precisionUpdateFactor;
         storm::RationalNumber maxVerificationIterationFactor;
-        bool relevantValuesForPrecisionUpdate;
         storm::RationalNumber upperBoundGuessingFactor;
         uint64_t upperBoundOnlyIterations;
         bool noTerminationGuaranteeMinimumMethod;

src/storm/modelchecker/CheckTask.h

@@ -207,6 +207,14 @@ namespace storm {
             bool isQualitativeSet() const {
                 return qualitative;
             }
+
+            /*!
+             * sets whether the computation only needs to be performed qualitatively, because the values will only
+             * be compared to 0/1.
+             */
+            void setQualitative(bool value) {
+                qualitative = value;
+            }
             
             /*!
              * Sets whether to produce schedulers (if supported).

src/storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessor.cpp

@@ -343,8 +343,6 @@ namespace storm {
     
                 template<typename SparseModelType>
                 void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessTimeOperatorFormula(storm::logic::TimeOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
-                    // Time formulas are only supported for Markov automata
-                    STORM_LOG_THROW(data.model->isOfType(storm::models::ModelType::MarkovAutomaton), storm::exceptions::InvalidPropertyException, "Time operator formulas are only supported for Markov automata.");
                     
                     data.objectives.back()->lowerResultBound = storm::utility::zero<ValueType>();
                     
@@ -474,11 +472,15 @@ namespace storm {
                                 }
                             }
                             data.model->addRewardModel(rewardModelName, std::move(objectiveRewards));
-                        } else if (formula.isReachabilityTimeFormula() && data.model->isOfType(storm::models::ModelType::MarkovAutomaton)) {
+                        } else if (formula.isReachabilityTimeFormula()) {
                             
-                            // build stateAction reward vector that only gives reward for relevant states
+                            // build state reward vector that only gives reward for relevant states
                             std::vector<typename SparseModelType::ValueType> timeRewards(data.model->getNumberOfStates(), storm::utility::zero<typename SparseModelType::ValueType>());
-                            storm::utility::vector::setVectorValues(timeRewards, dynamic_cast<storm::models::sparse::MarkovAutomaton<typename SparseModelType::ValueType> const&>(*data.model).getMarkovianStates() & reachableFromInit, storm::utility::one<typename SparseModelType::ValueType>());
+                            if (data.model->isOfType(storm::models::ModelType::MarkovAutomaton)) {
+                                storm::utility::vector::setVectorValues(timeRewards, dynamic_cast<storm::models::sparse::MarkovAutomaton<typename SparseModelType::ValueType> const&>(*data.model).getMarkovianStates() & reachableFromInit, storm::utility::one<typename SparseModelType::ValueType>());
+                            } else {
+                                storm::utility::vector::setVectorValues(timeRewards, reachableFromInit, storm::utility::one<typename SparseModelType::ValueType>());
+                            }
                             data.model->addRewardModel(rewardModelName, typename SparseModelType::RewardModelType(std::move(timeRewards)));
                         } else {
                             STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The formula " << formula << " neither considers reachability probabilities nor reachability rewards " << (data.model->isOfType(storm::models::ModelType::MarkovAutomaton) ?  "nor reachability time" : "") << ". This is not supported.");

src/storm/settings/modules/OviSolverSettings.cpp

@@ -14,7 +14,6 @@ namespace storm {
             const std::string OviSolverSettings::moduleName = "ovi";
             const std::string OviSolverSettings::precisionUpdateFactorOptionName = "precision-update-factor";
             const std::string OviSolverSettings::maxVerificationIterationFactorOptionName = "max-verification-iter-factor";
-            const std::string OviSolverSettings::useRelevantValuesForPrecisionUpdateOptionName = "use-relevant-values";
             const std::string OviSolverSettings::upperBoundGuessingFactorOptionName = "upper-bound-factor";
             const std::string OviSolverSettings::upperBoundOnlyIterationsOptionName = "check-upper-only-iter";
             const std::string OviSolverSettings::useNoTerminationGuaranteeMinimumMethodOptionName = "no-termination-guarantee";
@@ -23,9 +22,7 @@ namespace storm {
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, precisionUpdateFactorOptionName, false, "Sets with which factor the precision of the inner value iteration is updated.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("factor", "The factor.").setDefaultValueDouble(0.4).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
                 
-                this->addOption(storm::settings::OptionBuilder(moduleName, useRelevantValuesForPrecisionUpdateOptionName, false, "Sets whether the precision of the inner value iteration is only based on the relevant values (i.e. initial states).").setIsAdvanced().build());
-
-                this->addOption(storm::settings::OptionBuilder(moduleName, maxVerificationIterationFactorOptionName, false, "Controls how many verification iterations are performed before guessing a new upper bound.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("factor", "The factor.").setDefaultValueDouble(0.1).addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterValidator(0.0)).build()).build());
+                this->addOption(storm::settings::OptionBuilder(moduleName, maxVerificationIterationFactorOptionName, false, "Controls how many verification iterations are performed before guessing a new upper bound.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("factor", "The factor.").setDefaultValueDouble(1.0).addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterValidator(0.0)).build()).build());
 
                 this->addOption(storm::settings::OptionBuilder(moduleName, upperBoundGuessingFactorOptionName, false, "Sets with which factor the precision is multiplied to guess the upper bound.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("factor", "The factor.").setDefaultValueDouble(1.0).addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterValidator(0.0)).build()).build());
 
@@ -42,10 +39,6 @@ namespace storm {
                 return this->getOption(maxVerificationIterationFactorOptionName).getArgumentByName("factor").getValueAsDouble();
             }
             
-            bool OviSolverSettings::useRelevantValuesForPrecisionUpdate() const {
-                return this->getOption(useRelevantValuesForPrecisionUpdateOptionName).getHasOptionBeenSet();
-            }
-
             double OviSolverSettings::getUpperBoundGuessingFactor() const {
                 return this->getOption(upperBoundGuessingFactorOptionName).getArgumentByName("factor").getValueAsDouble();
             }

src/storm/settings/modules/OviSolverSettings.h

@@ -19,8 +19,6 @@ namespace storm {
                 
                 double getMaxVerificationIterationFactor() const;
                 
-                bool useRelevantValuesForPrecisionUpdate() const;
-
                 double getUpperBoundGuessingFactor() const;
 
                 uint64_t getUpperBoundOnlyIterations() const;
@@ -33,7 +31,6 @@ namespace storm {
             private:
                 static const std::string precisionUpdateFactorOptionName;
                 static const std::string maxVerificationIterationFactorOptionName;
-                static const std::string useRelevantValuesForPrecisionUpdateOptionName;
                 static const std::string upperBoundGuessingFactorOptionName;
                 static const std::string upperBoundOnlyIterationsOptionName;
                 static const std::string useNoTerminationGuaranteeMinimumMethodOptionName;

src/storm/solver/AbstractEquationSolver.cpp

@@ -66,6 +66,11 @@ namespace storm {
             return relevantValues.get();
         }
         
+        template<typename ValueType>
+        boost::optional<storm::storage::BitVector> const& AbstractEquationSolver<ValueType>::getOptionalRelevantValues() const {
+            return relevantValues;
+        }
+        
         template<typename ValueType>
         void AbstractEquationSolver<ValueType>::setRelevantValues(storm::storage::BitVector&& relevantValues) {
             this->relevantValues = std::move(relevantValues);

src/storm/solver/AbstractEquationSolver.h

@@ -52,6 +52,7 @@ namespace storm {
              * Retrieves the relevant values (if there are any).
              */
             storm::storage::BitVector const& getRelevantValues() const;
+            boost::optional<storm::storage::BitVector> const& getOptionalRelevantValues() const;
             
             /*!
              * Sets the relevant values.

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -2,7 +2,6 @@
 #include <limits>
 
 #include "storm/solver/IterativeMinMaxLinearEquationSolver.h"
-#include "storm/solver/helper/OptimisticValueIterationHelper.h"
 
 #include "storm/environment/solver/MinMaxSolverEnvironment.h"
 #include "storm/environment/solver/OviSolverEnvironment.h"
@@ -369,63 +368,28 @@ namespace storm {
                 return true;
             }
             
-            if (!this->multiplierA) {
-                this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
-            }
-
             if (!auxiliaryRowGroupVector) {
                 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
             }
-            if (!auxiliaryRowGroupVector2) {
-                auxiliaryRowGroupVector2 = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
+            if (!optimisticValueIterationHelper) {
+                optimisticValueIterationHelper = std::make_unique<storm::solver::helper::OptimisticValueIterationHelper<ValueType>>(*this->A);
             }
 
-            // By default, we can not provide any guarantee
-            SolverGuarantee guarantee = SolverGuarantee::None;
-            // Get handle to multiplier.
-            storm::solver::Multiplier<ValueType> const &multiplier = *this->multiplierA;
-            // Allow aliased multiplications.
-            storm::solver::MultiplicationStyle multiplicationStyle = env.solver().minMax().getMultiplicationStyle();
-            bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel;
-
-            boost::optional<storm::storage::BitVector> relevantValues;
-            if (this->hasRelevantValues()) {
-                relevantValues = this->getRelevantValues();
-            }
+            storm::solver::helper::OptimisticValueIterationHelper<ValueType> helper(*this->A);
             
             // x has to start with a lower bound.
             this->createLowerBoundsVector(x);
 
             std::vector<ValueType>* lowerX = &x;
             std::vector<ValueType>* upperX = auxiliaryRowGroupVector.get();
-            std::vector<ValueType>* auxVector = auxiliaryRowGroupVector2.get();
 
-            this->startMeasureProgress();
             
-            storm::solver::helper::OptimisticValueIterationHelper<ValueType> helper(*this->A);
-            auto statusIters = helper.solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, b,
-                    [&] (std::vector<ValueType>*& y, std::vector<ValueType>*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI) {
-                        this->showProgressIterative(i);
-                        auto result = performValueIteration(env, dir, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle);
-                        return std::make_pair(result.iterations, result.status);
-                    },
-                    [&] (std::vector<ValueType>* y, std::vector<ValueType>* yPrime, uint64_t const& i) {
-                        this->showProgressIterative(i);
-                        if (useGaussSeidelMultiplication) {
-                            // Copy over the current vectors so we can modify them in-place.
-                            // This is necessary as we want to compare the new values with the current ones.
-                            *yPrime = *y;
-                            multiplier.multiplyAndReduceGaussSeidel(env, dir, *y, &b);
-                        } else {
-                            multiplier.multiplyAndReduce(env, dir, *y, &b, *yPrime);
-                            std::swap(y, yPrime);
-                        }
-                    },
-                    env.solver().minMax().getRelativeTerminationCriterion(),
-                    storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()),
-                    env.solver().minMax().getMaximalNumberOfIterations(),
-                    dir,
-                    relevantValues);
+            auto statusIters = helper.solveEquations(env, lowerX, upperX, b,
+                                                     env.solver().minMax().getRelativeTerminationCriterion(),
+                                                     storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()),
+                                                     env.solver().minMax().getMaximalNumberOfIterations(),
+                                                     dir,
+                                                     this->getOptionalRelevantValues());
             auto two = storm::utility::convertNumber<ValueType>(2.0);
             storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(*lowerX, *upperX, x, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; });
 
@@ -434,8 +398,7 @@ namespace storm {
             // If requested, we store the scheduler for retrieval.
             if (this->isTrackSchedulerSet()) {
                 this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
-                this->multiplierA->multiplyAndReduce(env, dir, x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get());
-                this->multiplierA->multiplyAndReduce(env, dir, x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get());
+                this->A->multiplyAndReduce(dir, this->A->getRowGroupIndices(), x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get());
             }
 
             if (!this->isCachingEnabled()) {
@@ -1109,6 +1072,7 @@ namespace storm {
             auxiliaryRowGroupVector.reset();
             auxiliaryRowGroupVector2.reset();
             soundValueIterationHelper.reset();
+            optimisticValueIterationHelper.reset();
             StandardMinMaxLinearEquationSolver<ValueType>::clearCache();
         }
         

src/storm/solver/IterativeMinMaxLinearEquationSolver.h

@@ -8,6 +8,7 @@
 #include "storm/solver/Multiplier.h"
 #include "storm/solver/StandardMinMaxLinearEquationSolver.h"
 #include "storm/solver/helper/SoundValueIterationHelper.h"
+#include "storm/solver/helper/OptimisticValueIterationHelper.h"
 
 #include "storm/solver/SolverStatus.h"
 
@@ -85,6 +86,7 @@ namespace storm {
             mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRowGroupVector; // A.rowGroupCount() entries
             mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRowGroupVector2; // A.rowGroupCount() entries
             mutable std::unique_ptr<storm::solver::helper::SoundValueIterationHelper<ValueType>> soundValueIterationHelper;
+            mutable std::unique_ptr<storm::solver::helper::OptimisticValueIterationHelper<ValueType>> optimisticValueIterationHelper;
             
         };
         

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -643,62 +643,26 @@ namespace storm {
                 return true;
             }
             
-            if (!this->multiplier) {
-                this->multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
-            }
-
             if (!this->cachedRowVector) {
                 this->cachedRowVector = std::make_unique<std::vector<ValueType>>(this->A->getRowCount());
             }
-            if (!this->cachedRowVector2) {
-                this->cachedRowVector2 = std::make_unique<std::vector<ValueType>>(this->A->getRowCount());
+            if (!optimisticValueIterationHelper) {
+                optimisticValueIterationHelper = std::make_unique<storm::solver::helper::OptimisticValueIterationHelper<ValueType>>(*this->A);
             }
 
-            // By default, we can not provide any guarantee
-            SolverGuarantee guarantee = SolverGuarantee::None;
-            // Get handle to multiplier.
-            storm::solver::Multiplier<ValueType> const &multiplier = *this->multiplier;
-            // Allow aliased multiplications.
-            storm::solver::MultiplicationStyle multiplicationStyle = env.solver().native().getPowerMethodMultiplicationStyle();
-            bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel;
-
-            boost::optional<storm::storage::BitVector> relevantValues;
-            if (this->hasRelevantValues()) {
-                relevantValues = this->getRelevantValues();
-            }
-            
             // x has to start with a lower bound.
             this->createLowerBoundsVector(x);
 
             std::vector<ValueType>* lowerX = &x;
             std::vector<ValueType>* upperX = this->cachedRowVector.get();
-            std::vector<ValueType>* auxVector = this->cachedRowVector2.get();
 
-            this->startMeasureProgress();
             storm::solver::helper::OptimisticValueIterationHelper<ValueType> helper(*this->A);
-            auto statusIters = helper.solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, b,
-                    [&] (std::vector<ValueType>*& y, std::vector<ValueType>*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI) {
-                        this->showProgressIterative(i);
-                        auto result = performPowerIteration(env, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle);
-                        return std::make_pair(result.iterations, result.status);
-                    },
-                    [&] (std::vector<ValueType>* y, std::vector<ValueType>* yPrime, uint64_t const& i) {
-                        this->showProgressIterative(i);
-                        if (useGaussSeidelMultiplication) {
-                            // Copy over the current vectors so we can modify them in-place.
-                            // This is necessary as we want to compare the new values with the current ones.
-                            *yPrime = *y;
-                            multiplier.multiplyGaussSeidel(env, *y, &b);
-                        } else {
-                            multiplier.multiply(env, *y, &b, *yPrime);
-                            std::swap(y, yPrime);
-                        }
-                    },
-                    env.solver().native().getRelativeTerminationCriterion(),
-                    storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision()),
-                    env.solver().native().getMaximalNumberOfIterations(),
-                    boost::none, // No optimization dir
-                    relevantValues);
+            auto statusIters = helper.solveEquations(env, lowerX, upperX, b,
+                                                     env.solver().native().getRelativeTerminationCriterion(),
+                                                     storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision()),
+                                                     env.solver().native().getMaximalNumberOfIterations(),
+                                                     boost::none, // No optimization dir
+                                                     this->getOptionalRelevantValues());
             auto two = storm::utility::convertNumber<ValueType>(2.0);
             storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(*lowerX, *upperX, x, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; });
             this->logIterations(statusIters.first == SolverStatus::Converged, statusIters.first == SolverStatus::TerminatedEarly, statusIters.second);
@@ -1066,6 +1030,7 @@ namespace storm {
             walkerChaeData.reset();
             multiplier.reset();
             soundValueIterationHelper.reset();
+            optimisticValueIterationHelper.reset();
             LinearEquationSolver<ValueType>::clearCache();
         }
         

src/storm/solver/NativeLinearEquationSolver.h

@@ -9,6 +9,7 @@
 #include "storm/solver/NativeMultiplier.h"
 #include "storm/solver/SolverStatus.h"
 #include "storm/solver/helper/SoundValueIterationHelper.h"
+#include "storm/solver/helper/OptimisticValueIterationHelper.h"
 
 #include "storm/utility/NumberTraits.h"
 
@@ -96,6 +97,7 @@ namespace storm {
             // cached auxiliary data
             mutable std::unique_ptr<std::vector<ValueType>> cachedRowVector2; // A.getRowCount() rows
             mutable std::unique_ptr<storm::solver::helper::SoundValueIterationHelper<ValueType>> soundValueIterationHelper;
+            mutable std::unique_ptr<storm::solver::helper::OptimisticValueIterationHelper<ValueType>> optimisticValueIterationHelper;
             
             struct JacobiDecomposition {
                 JacobiDecomposition(Environment const& env, storm::storage::SparseMatrix<ValueType> const& A);

src/storm/solver/TopologicalLinearEquationSolver.cpp

@@ -144,11 +144,12 @@ namespace storm {
             }
             
             if (hasDiagonalEntry) {
-                    if (storm::utility::isZero(denominator)) {
-                        STORM_LOG_THROW(storm::utility::isZero(xi), storm::exceptions::InvalidOperationException, "The equation system has no solution.");
-                    } else {
-                        xi /= denominator;
-                    }
+                STORM_LOG_WARN_COND_DEBUG(storm::NumberTraits<ValueType>::IsExact || !storm::utility::isAlmostZero(denominator) || storm::utility::isZero(denominator), "State " << sccState << " has a selfloop with probability '1-(" << denominator << ")'. This could be an indication for numerical issues.");
+                if (storm::utility::isZero(denominator)) {
+                    STORM_LOG_THROW(storm::utility::isZero(xi), storm::exceptions::InvalidOperationException, "The equation system has no solution.");
+                } else {
+                    xi /= denominator;
+                }
             }
             return true;
         }

src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp

@@ -155,6 +155,7 @@ namespace storm {
                     }
                 }
                 if (hasDiagonalEntry) {
+                    STORM_LOG_WARN_COND_DEBUG(storm::NumberTraits<ValueType>::IsExact || !storm::utility::isAlmostZero(denominator) || storm::utility::isZero(denominator), "State " << sccState << " has a selfloop with probability '1-(" << denominator << ")'. This could be an indication for numerical issues.");
                     if (storm::utility::isZero(denominator)) {
                         // In this case we have a selfloop on this state. This can never an optimal choice:
                         // When minimizing, we are looking for the largest fixpoint (which will never be attained by this action)

src/storm/solver/helper/OptimisticValueIterationHelper.cpp

@@ -15,56 +15,9 @@ namespace storm {
             namespace oviinternal {
     
                 template<typename ValueType>
-                ValueType computeMaxAbsDiff(std::vector<ValueType> const& allOldValues, std::vector<ValueType> const& allNewValues, storm::storage::BitVector const& relevantValues) {
-                    ValueType result = storm::utility::zero<ValueType>();
-                    for (auto value : relevantValues) {
-                        result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allNewValues[value] - allOldValues[value]));
-                    }
-                    return result;
-                }
-        
-                template<typename ValueType>
-                ValueType computeMaxAbsDiff(std::vector<ValueType> const& allOldValues, std::vector<ValueType> const& allNewValues) {
-                    ValueType result = storm::utility::zero<ValueType>();
-                    for (uint64_t i = 0; i < allOldValues.size(); ++i) {
-                        result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allNewValues[i] - allOldValues[i]));
-                    }
-                    return result;
-                }
-            
-                template<typename ValueType>
-                ValueType computeMaxRelDiff(std::vector<ValueType> const& allOldValues, std::vector<ValueType> const& allNewValues, storm::storage::BitVector const& relevantValues) {
-                    ValueType result = storm::utility::zero<ValueType>();
-                    for (auto const& i : relevantValues) {
-                        STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector.");
-                        if (!storm::utility::isZero(allNewValues[i])) {
-                            result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allNewValues[i] - allOldValues[i]) / allNewValues[i]);
-                        }
-                    }
-                    return result;
-                }
-                
-                template<typename ValueType>
-                ValueType computeMaxRelDiff(std::vector<ValueType> const& allOldValues, std::vector<ValueType> const& allNewValues) {
-                    ValueType result = storm::utility::zero<ValueType>();
-                    for (uint64_t i = 0; i < allOldValues.size(); ++i) {
-                        STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector.");
-                        if (!storm::utility::isZero(allNewValues[i])) {
-                            result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allNewValues[i] - allOldValues[i]) / allNewValues[i]);
-                        }
-                    }
-                    return result;
-                }
-    
-                template<typename ValueType>
-                ValueType updateIterationPrecision(storm::Environment const& env, std::vector<ValueType> const& currentX, std::vector<ValueType> const& newX, bool const& relative, boost::optional<storm::storage::BitVector> const& relevantValues) {
+                ValueType updateIterationPrecision(storm::Environment const& env, ValueType const& diff) {
                     auto factor = storm::utility::convertNumber<ValueType>(env.solver().ovi().getPrecisionUpdateFactor());
-                    bool useRelevant = relevantValues.is_initialized() && env.solver().ovi().useRelevantValuesForPrecisionUpdate();
-                    if (relative) {
-                        return (useRelevant ? computeMaxRelDiff(newX, currentX, relevantValues.get()) : computeMaxRelDiff(newX, currentX)) * factor;
-                    } else {
-                        return (useRelevant ? computeMaxAbsDiff(newX, currentX, relevantValues.get()) : computeMaxAbsDiff(newX, currentX)) * factor;
-                    }
+                    return factor * diff;
                 }
     
                 template<typename ValueType>
@@ -78,7 +31,7 @@ namespace storm {
                 }
                 
                 template <typename ValueType>
-                UpperBoundIterator<ValueType>::UpperBoundIterator(storm::storage::SparseMatrix<ValueType> const& matrix) {
+                IterationHelper<ValueType>::IterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix) {
                     STORM_LOG_THROW(static_cast<uint64_t>(std::numeric_limits<uint32_t>::max()) > matrix.getRowCount() + 1, storm::exceptions::NotSupportedException, "Matrix dimensions too large.");
                     STORM_LOG_THROW(static_cast<uint64_t>(std::numeric_limits<uint32_t>::max()) > matrix.getEntryCount(), storm::exceptions::NotSupportedException, "Matrix dimensions too large.");
                     matrixValues.reserve(matrix.getNonzeroEntryCount());
@@ -98,22 +51,118 @@ namespace storm {
                 }
                 
                 template <typename ValueType>
-                typename UpperBoundIterator<ValueType>::IterateResult UpperBoundIterator<ValueType>::iterate(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
-                    return iterateInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew);
+                ValueType IterationHelper<ValueType>::singleIterationWithDiff(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) {
+                    return singleIterationWithDiffInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, computeRelativeDiff);
+                }
+                
+                template <typename ValueType>
+                ValueType IterationHelper<ValueType>::singleIterationWithDiff(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) {
+                    if (minimize(dir)) {
+                        return singleIterationWithDiffInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, computeRelativeDiff);
+                    } else {
+                        return singleIterationWithDiffInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, computeRelativeDiff);
+                    }
+                }
+                
+                template <typename ValueType>
+                template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
+                ValueType IterationHelper<ValueType>::singleIterationWithDiffInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) {
+                    STORM_LOG_ASSERT(x.size() > 0, "Empty equation system not expected.");
+                    ValueType diff = storm::utility::zero<ValueType>();
+                    
+                    IndexType i = x.size();
+                    while (i > 0) {
+                        --i;
+                        ValueType newXi = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x);
+                        ValueType& oldXi = x[i];
+                        if (computeRelativeDiff) {
+                            if (storm::utility::isZero(newXi)) {
+                                if (storm::utility::isZero(oldXi)) {
+                                    // this operation has no effect:
+                                    // diff = std::max(diff, storm::utility::zero<ValueType>());
+                                } else {
+                                    diff = std::max(diff, storm::utility::one<ValueType>());
+                                }
+                            } else {
+                                diff = std::max(diff, storm::utility::abs<ValueType>((newXi - oldXi) / newXi));
+                            }
+                        } else {
+                            diff = std::max(diff, storm::utility::abs<ValueType>(newXi - oldXi));
+                        }
+                        oldXi = std::move(newXi);
+                    }
+                    return diff;
                 }
                 
                 template <typename ValueType>
-                typename UpperBoundIterator<ValueType>::IterateResult UpperBoundIterator<ValueType>::iterate(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
+                uint64_t IterationHelper<ValueType>::repeatedIterate(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative) {
                     if (minimize(dir)) {
-                        return iterateInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew);
+                        return repeatedIterateInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, precision, relative);
                     } else {
-                        return iterateInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, takeMinOfOldAndNew);
+                        return repeatedIterateInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, precision, relative);
                     }
                 }
                 
+                template <typename ValueType>
+                uint64_t IterationHelper<ValueType>::repeatedIterate(std::vector<ValueType>& x, const std::vector<ValueType>& b, ValueType precision, bool relative) {
+                    return repeatedIterateInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, precision, relative);
+                }
+                
                 template <typename ValueType>
                 template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
-                typename UpperBoundIterator<ValueType>::IterateResult UpperBoundIterator<ValueType>::iterateInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
+                uint64_t IterationHelper<ValueType>::repeatedIterateInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative) {
+                    // Do a backwards gauss-seidel style iteration
+                    bool convergence = true;
+                    IndexType i = x.size();
+                    while (i > 0) {
+                        --i;
+                        ValueType newXi = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x);
+                        ValueType& oldXi = x[i];
+                        // Check if we converged
+                        if (relative) {
+                            if (storm::utility::isZero(oldXi)) {
+                                if (!storm::utility::isZero(newXi)) {
+                                    convergence = false;
+                                    break;
+                                }
+                            } else if (storm::utility::abs<ValueType>((newXi - oldXi) / oldXi) > precision) {
+                                convergence = false;
+                                break;
+                            }
+                        } else {
+                            if (storm::utility::abs<ValueType>((newXi - oldXi)) > precision) {
+                                convergence = false;
+                                break;
+                            }
+                        }
+                    }
+                    if (!convergence) {
+                        // we now know that we did not converge. We still need to set the remaining values
+                        while (i > 0) {
+                            --i;
+                            x[i] = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x);
+                        }
+                    }
+                    return convergence;
+                }
+                
+                template <typename ValueType>
+                typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpper(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
+                    if (minimize(dir)) {
+                        return iterateUpperInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew);
+                    } else {
+                        return iterateUpperInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, takeMinOfOldAndNew);
+                    }
+                }
+                
+                template <typename ValueType>
+                typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpper(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
+                    return iterateUpperInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew);
+                }
+                
+                template <typename ValueType>
+                template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
+                typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpperInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) {
                     // For each row compare the new upper bound candidate with the old one
                     bool newUpperBoundAlwaysHigherEqual = true;
                     bool newUpperBoundAlwaysLowerEqual = true;
@@ -150,7 +199,7 @@ namespace storm {
                 }
                 
                 template <typename ValueType>
-                ValueType UpperBoundIterator<ValueType>::multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector<ValueType> const& x) {
+                ValueType IterationHelper<ValueType>::multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector<ValueType> const& x) {
                     assert(rowIndex < rowIndications.size());
                     ValueType xRes = bi;
                     
@@ -165,7 +214,7 @@ namespace storm {
                 
                 template <typename ValueType>
                 template<storm::solver::OptimizationDirection Dir>
-                ValueType UpperBoundIterator<ValueType>::multiplyRowGroup(IndexType const& rowGroupIndex, std::vector<ValueType> const& b, std::vector<ValueType> const& x) {
+                ValueType IterationHelper<ValueType>::multiplyRowGroup(IndexType const& rowGroupIndex, std::vector<ValueType> const& b, std::vector<ValueType> const& x) {
                     STORM_LOG_ASSERT(rowGroupIndices != nullptr, "No row group indices available.");
                     auto row = (*rowGroupIndices)[rowGroupIndex];
                     auto const& groupEnd = (*rowGroupIndices)[rowGroupIndex + 1];
@@ -180,24 +229,21 @@ namespace storm {
             }
             
             template<typename ValueType>
-            OptimisticValueIterationHelper<ValueType>::OptimisticValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix) : upperBoundIterator(matrix) {
+            OptimisticValueIterationHelper<ValueType>::OptimisticValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix) : iterationHelper(matrix) {
                     // Intentionally left empty.
-                }
+            }
                 
             template<typename ValueType>
-                std::pair<SolverStatus, uint64_t> OptimisticValueIterationHelper<ValueType>::solveEquationsOptimisticValueIteration(Environment const& env, std::vector<ValueType>* lowerX, std::vector<ValueType>* upperX, std::vector<ValueType>* auxVector, std::vector<ValueType> const& b, ValueIterationCallBackType const& valueIterationCallback, SingleIterationCallBackType const& singleIterationCallback, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional<storm::solver::OptimizationDirection> dir, boost::optional<storm::storage::BitVector> relevantValues) {
+                std::pair<SolverStatus, uint64_t> OptimisticValueIterationHelper<ValueType>::solveEquations(Environment const& env, std::vector<ValueType>* lowerX, std::vector<ValueType>* upperX, std::vector<ValueType> const& b, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional<storm::solver::OptimizationDirection> dir, boost::optional<storm::storage::BitVector> const& relevantValues) {
                 STORM_LOG_ASSERT(lowerX->size() == upperX->size(), "Dimension missmatch.");
-                STORM_LOG_ASSERT(lowerX->size() == auxVector->size(), "Dimension missmatch.");
                 
                 // As we will shuffle pointers around, let's store the original positions here.
                 std::vector<ValueType>* initLowerX = lowerX;
                 std::vector<ValueType>* initUpperX = upperX;
-                std::vector<ValueType>* initAux = auxVector;
                 
                 uint64_t overallIterations = 0;
                 uint64_t lastValueIterationIterations = 0;
                 uint64_t currentVerificationIterations = 0;
-                uint64_t valueIterationInvocations = 0;
                 
                 // Get some parameters for the algorithm
                 // 2
@@ -215,135 +261,101 @@ namespace storm {
                 SolverStatus status = SolverStatus::InProgress;
     
                 while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) {
-
                     // Perform value iteration until convergence
-                    ++valueIterationInvocations;
-                    auto result = valueIterationCallback(lowerX, auxVector, iterationPrecision, relative, overallIterations, maxOverallIterations);
-                    lastValueIterationIterations = result.first;
-                    overallIterations += result.first;
+                    lastValueIterationIterations = dir ? iterationHelper.repeatedIterate(dir.get(), *lowerX, b, iterationPrecision, relative) : iterationHelper.repeatedIterate(*lowerX, b, iterationPrecision, relative);
+                    overallIterations += lastValueIterationIterations;
 
-                    if (result.second != SolverStatus::Converged) {
-                        status = result.second;
-                    } else {
-                        bool intervalIterationNeeded = false;
-                        currentVerificationIterations = 0;
+                    bool intervalIterationNeeded = false;
+                    currentVerificationIterations = 0;
 
-                        if (relative) {
-                            oviinternal::guessUpperBoundRelative(*lowerX, *upperX, relativeBoundGuessingScaler);
-                        } else {
-                            oviinternal::guessUpperBoundAbsolute(*lowerX, *upperX, precision);
-                        }
+                    if (relative) {
+                        oviinternal::guessUpperBoundRelative(*lowerX, *upperX, relativeBoundGuessingScaler);
+                    } else {
+                        oviinternal::guessUpperBoundAbsolute(*lowerX, *upperX, precision);
+                    }
 
-                        bool cancelGuess = false;
-                        while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) {
-                            ++overallIterations;
-                            ++currentVerificationIterations;
-                            // Perform value iteration stepwise for lower bound and guessed upper bound
+                    bool cancelGuess = false;
+                    while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) {
+                        ++overallIterations;
+                        ++currentVerificationIterations;
+                        // Perform value iteration stepwise for lower bound and guessed upper bound
 
-                            // Upper bound iteration
-                            auto upperBoundIterResult = dir ? upperBoundIterator.iterate(dir.get(), *upperX, b, !noTerminationGuarantee) : upperBoundIterator.iterate(*upperX, b, !noTerminationGuarantee);
+                        // Upper bound iteration
+                        auto upperBoundIterResult = dir ? iterationHelper.iterateUpper(dir.get(), *upperX, b, !noTerminationGuarantee) : iterationHelper.iterateUpper(*upperX, b, !noTerminationGuarantee);
 
-                            if (upperBoundIterResult == oviinternal::UpperBoundIterator<ValueType>::IterateResult::AlwaysHigherOrEqual) {
-                                // All values moved up (and did not stay the same)
-                                // That means the guess for an upper bound is actually a lower bound
-                                iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *upperX, relative, relevantValues);
-                                // We assume to have a single fixed point. We can thus safely set the new lower bound, to the wrongly guessed upper bound
-                                // Set lowerX to the upper bound candidate
-                                std::swap(lowerX, upperX);
-                                break;
-                            } else if (upperBoundIterResult == oviinternal::UpperBoundIterator<ValueType>::IterateResult::AlwaysLowerOrEqual) {
-                                // All values moved down (and stayed not the same)
-                                // This is a valid upper bound. We still need to check the precision.
-                                // We can safely use twice the requested precision, as we calculate the center of both vectors
-                                bool reachedPrecision;
-                                if (relevantValues) {
-                                    reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, relevantValues.get(), doublePrecision, relative);
-                                } else {
-                                    reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, doublePrecision, relative);
-                                }
-                                if (reachedPrecision) {
-                                    status = SolverStatus::Converged;
-                                    break;
-                                } else {
-                                    // From now on, we keep updating both bounds
-                                    intervalIterationNeeded = true;
-                                }
-                            // The following case below covers that both vectors (old and new) are equal.
-                            // Theoretically, this means that the precise fixpoint has been reached. However, numerical instabilities can be tricky and this detection might be incorrect (see the haddad-monmege model).
-                            // We therefore disable it. It is very unlikely that we guessed the right fixpoint anyway.
-                            //} else if (upperBoundIterResult == oviinternal::UpperBoundIterator<ValueType>::IterateResult::Equal) {
-                                // In this case, the guessed upper bound is the precise fixpoint
-                            //    status = SolverStatus::Converged;
-                            //    std::swap(lowerX, auxVector);
-                            //    break;
+                        if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::AlwaysHigherOrEqual) {
+                            // All values moved up (and did not stay the same)
+                            // That means the guess for an upper bound is actually a lower bound
+                            auto diff = dir ? iterationHelper.singleIterationWithDiff(dir.get(), *upperX, b, relative) : iterationHelper.singleIterationWithDiff(*upperX, b, relative);
+                            iterationPrecision = oviinternal::updateIterationPrecision(env, diff);
+                            // We assume to have a single fixed point. We can thus safely set the new lower bound, to the wrongly guessed upper bound
+                            // Set lowerX to the upper bound candidate
+                            std::swap(lowerX, upperX);
+                            break;
+                        } else if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::AlwaysLowerOrEqual) {
+                            // All values moved down (and stayed not the same)
+                            // This is a valid upper bound. We still need to check the precision.
+                            // We can safely use twice the requested precision, as we calculate the center of both vectors
+                            bool reachedPrecision;
+                            if (relevantValues) {
+                                reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, relevantValues.get(), doublePrecision, relative);
+                            } else {
+                                reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, doublePrecision, relative);
                             }
-
-                            // Check whether we tried this guess for too long
-                            ValueType scaledIterationCount = storm::utility::convertNumber<ValueType>(currentVerificationIterations) * storm::utility::convertNumber<ValueType>(env.solver().ovi().getMaxVerificationIterationFactor());
-                            if (!intervalIterationNeeded && scaledIterationCount >= storm::utility::convertNumber<ValueType>(lastValueIterationIterations)) {
-                                cancelGuess = true;
-                                // In this case we will make one more iteration on the lower bound (mainly to obtain a new iterationPrecision)
+                            if (reachedPrecision) {
+                                status = SolverStatus::Converged;
+                                break;
+                            } else {
+                                // From now on, we keep updating both bounds
+                                intervalIterationNeeded = true;
                             }
+                        // The following case below covers that both vectors (old and new) are equal.
+                        // Theoretically, this means that the precise fixpoint has been reached. However, numerical instabilities can be tricky and this detection might be incorrect (see the haddad-monmege model).
+                        // We therefore disable it. It is very unlikely that we guessed the right fixpoint anyway.
+                        //} else if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::Equal) {
+                            // In this case, the guessed upper bound is the precise fixpoint
+                        //    status = SolverStatus::Converged;
+                        //    break;
+                        }
 
-                            // Lower bound iteration (only if needed)
-                            if (cancelGuess || intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) {
-                                singleIterationCallback(lowerX, auxVector, overallIterations);
-                                // At this point, auxVector contains the old values for the lower bound whereas lowerX contains the new ones.
+                        // Check whether we tried this guess for too long
+                        ValueType scaledIterationCount = storm::utility::convertNumber<ValueType>(currentVerificationIterations) * storm::utility::convertNumber<ValueType>(env.solver().ovi().getMaxVerificationIterationFactor());
+                        if (!intervalIterationNeeded && scaledIterationCount * iterationPrecision >= storm::utility::one<ValueType>()) {
+                            cancelGuess = true;
+                            // In this case we will make one more iteration on the lower bound (mainly to obtain a new iterationPrecision)
+                        }
 
-                                // Check whether the upper and lower bounds have crossed, i.e., the upper bound is smaller than the lower bound.
-                                bool valuesCrossed = false;
-                                for (uint64_t i = 0; i < lowerX->size(); ++i) {
-                                    if ((*upperX)[i] < (*lowerX)[i]) {
-                                        valuesCrossed = true;
-                                        break;
-                                    }
-                                }
+                        // Lower bound iteration (only if needed)
+                        if (cancelGuess || intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) {
+                            auto diff = dir ? iterationHelper.singleIterationWithDiff(dir.get(), *lowerX, b, relative) : iterationHelper.singleIterationWithDiff(*lowerX, b, relative);
 
-                                if (cancelGuess || valuesCrossed) {
-                                    // A new guess is needed.
-                                    iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *lowerX, relative, relevantValues);
+                            // Check whether the upper and lower bounds have crossed, i.e., the upper bound is smaller than the lower bound.
+                            bool valuesCrossed = false;
+                            for (uint64_t i = 0; i < lowerX->size(); ++i) {
+                                if ((*upperX)[i] < (*lowerX)[i]) {
+                                    valuesCrossed = true;
                                     break;
                                 }
                             }
-                        }
-                        if (storm::utility::resources::isTerminate()) {
-                            status = SolverStatus::Aborted;
+
+                            if (cancelGuess || valuesCrossed) {
+                                // A new guess is needed.
+                                iterationPrecision = oviinternal::updateIterationPrecision(env, diff);
+                                break;
+                            }
                         }
                     }
+                    if (storm::utility::resources::isTerminate()) {
+                        status = SolverStatus::Aborted;
+                    }
                 } // end while
                 
-                // Swap the results into the output vectors.
-                if (initLowerX == lowerX) {
-                    // lowerX is already at the correct position. We still have to care for upperX
-                    if (initUpperX != upperX) {
-                        // UpperX is not at the correct position. It has to be at the auxVector
-                        assert(initAux == upperX);
-                        std::swap(*initUpperX, *initAux);
-                    }
-                } else if (initUpperX == upperX) {
-                    // UpperX is already at the correct position.
-                    // We already know that lowerX is at the wrong position. It has to be at the auxVector
-                    assert(initAux == lowerX);
-                    std::swap(*initLowerX, *initAux);
-                } else if (initAux == auxVector) {
-                    // We know that upperX and lowerX are swapped.
-                    assert(initLowerX == upperX);
+                // Swap the results into the output vectors (if necessary).
+                assert(initLowerX != lowerX || (initLowerX == lowerX && initUpperX == upperX));
+                if (initLowerX != lowerX) {
                     assert(initUpperX == lowerX);
-                    std::swap(*initUpperX, *initLowerX);
-                } else {
-                    // Now we know that all vectors are at the wrong position. There are only two possibilities left
-                    if (initLowerX == upperX) {
-                        assert(initUpperX == auxVector);
-                        assert(initAux == lowerX);
-                        std::swap(*initLowerX, *initAux);
-                        std::swap(*initUpperX, *initAux);
-                    } else {
-                        assert(initLowerX == auxVector);
-                        assert(initUpperX == lowerX);
-                        assert (initAux == upperX);
-                        std::swap(*initUpperX, *initAux);
-                        std::swap(*initLowerX, *initAux);
-                    }
+                    assert(initLowerX == upperX);
+                    lowerX->swap(*upperX);
                 }
                 
                 if (overallIterations > maxOverallIterations) {

src/storm/solver/helper/OptimisticValueIterationHelper.h

@@ -18,10 +18,10 @@ namespace storm {
             namespace oviinternal {
     
                 template <typename ValueType>
-                class UpperBoundIterator {
+                class IterationHelper {
                 public:
                     typedef uint32_t IndexType;
-                    UpperBoundIterator(storm::storage::SparseMatrix<ValueType> const& matrix);
+                    IterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix);
                     enum class IterateResult {
                         AlwaysHigherOrEqual,
                         AlwaysLowerOrEqual,
@@ -29,13 +29,26 @@ namespace storm {
                         Incomparable
                     };
                     
-                    IterateResult iterate(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
-                    IterateResult iterate(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
+                    /// performs a single iteration and returns the maximal difference between the iterations as well as the index where this difference happened
+                    ValueType singleIterationWithDiff(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff);
+                    ValueType singleIterationWithDiff(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff);
+                    
+                    /// returns the number of performed iterations
+                    uint64_t repeatedIterate(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative);
+                    uint64_t repeatedIterate(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative);
+                    
+                    /// Performs a single iteration for the upper bound
+                    IterateResult iterateUpper(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
+                    IterateResult iterateUpper(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
                 
                 private:
                     
                     template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
-                    IterateResult iterateInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
+                    ValueType singleIterationWithDiffInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff);
+                    template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
+                    uint64_t repeatedIterateInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative);
+                    template<bool HasRowGroups, storm::solver::OptimizationDirection Dir>
+                    IterateResult iterateUpperInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew);
                     ValueType multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector<ValueType> const& x);
                     template<storm::solver::OptimizationDirection Dir>
                     ValueType multiplyRowGroup(IndexType const& rowGroupIndex, std::vector<ValueType> const& b, std::vector<ValueType> const& x);
@@ -56,49 +69,20 @@ namespace storm {
             public:
                 OptimisticValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix);
     
-                /*!
-                 * Return type of value iteration callback.
-                 * The first component shall be the number of performed iterations, the second component is the final convergence status.
-                 */
-                typedef std::pair<uint64_t, storm::solver::SolverStatus> ValueIterationReturnType;
-                
-                /*!
-                 * Value iteration callback. Performs conventional value iteration for the given input.
-                 * @pre y points to a vector that contains starting values
-                 * @post y points to a vector that contains resulting values
-                 * @param yPrime points to auxiliary storage
-                 * @param precision is the target precision
-                 * @param relative sets whether relative precision is considered
-                 * @param maxI the maximal number of iterations to perform
-                 */
-                typedef std::function<ValueIterationReturnType(std::vector<ValueType>*& y, std::vector<ValueType>*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI)> ValueIterationCallBackType;
-                
-                /*!
-                 * Should perform a single value iteration step (using conventional value iteration).
-                 * @pre y points to a vector that contains starting values
-                 * @post y points to a vector that contains resulting values
-                 * @param yPrime points to auxiliary storage
-                 * @param currI the current iteration (can be used to display progress within the callback)
-                 */
-                typedef std::function<void(std::vector<ValueType>* y, std::vector<ValueType>* yPrime, uint64_t const& currI)> SingleIterationCallBackType;
-                
                 /*!
                  * @param env
                  * @param lowerX Needs to be some arbitrary lower bound on the actual values initially
                  * @param upperX Does not need to be an upper bound initially
-                 * @param auxVector auxiliary storage (same size as lowerX and upperX)
                  * @param b the values added to each matrix row (the b in A*x+b)
-                 * @param valueIterationCallback  Function that should perform standard value iteration on the input vector
-                 * @param singleIterationCallback Function that should perform a single value iteration step on the input vector e.g. ( x' = min/max(A*x + b))
                  * @param dir The optimization direction
                  * @param relevantValues If given, we only check the precision at the states with the given indices.
                  * @return The status upon termination as well as the number of iterations Also, the maximum (relative/absolute) difference between lowerX and upperX will be 2*epsilon
                  * with the provided precision parameters.
                  */
-                std::pair<SolverStatus, uint64_t> solveEquationsOptimisticValueIteration(Environment const& env, std::vector<ValueType>* lowerX, std::vector<ValueType>* upperX, std::vector<ValueType>* auxVector, std::vector<ValueType> const& b, ValueIterationCallBackType const& valueIterationCallback, SingleIterationCallBackType const& singleIterationCallback, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional<storm::solver::OptimizationDirection> dir, boost::optional<storm::storage::BitVector> relevantValues = boost::none);
+                std::pair<SolverStatus, uint64_t> solveEquations(Environment const& env, std::vector<ValueType>* lowerX, std::vector<ValueType>* upperX, std::vector<ValueType> const& b, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional<storm::solver::OptimizationDirection> dir, boost::optional<storm::storage::BitVector> const& relevantValues);
 
             private:
-                oviinternal::UpperBoundIterator<ValueType> upperBoundIterator;
+                oviinternal::IterationHelper<ValueType> iterationHelper;
                 
             };
         }

src/storm/storage/jani/EdgeContainer.cpp

@@ -67,9 +67,13 @@ namespace storm {
                     e.setTemplateEdge(map[e.getTemplateEdge()]);
                 }
             }
-
-
-
+        }
+        
+        EdgeContainer& EdgeContainer::operator=(EdgeContainer const& other) {
+            EdgeContainer otherCpy(other);
+            this->templates = std::move(otherCpy.templates);
+            this->edges = std::move(otherCpy.edges);
+            return *this;
         }
 
         void EdgeContainer::finalize(Model const& containingModel) {

src/storm/storage/jani/EdgeContainer.h

@@ -91,6 +91,7 @@ namespace storm {
 
             EdgeContainer() = default;
             EdgeContainer(EdgeContainer const& other);
+            EdgeContainer& operator=(EdgeContainer const& other);
 
             void clearConcreteEdges();
             std::vector<Edge> const& getConcreteEdges() const;

src/storm/storage/jani/TemplateEdgeContainer.cpp

@@ -3,10 +3,18 @@
 
 namespace storm {
     namespace jani {
-        TemplateEdgeContainer::TemplateEdgeContainer(TemplateEdgeContainer const &other) {
+        TemplateEdgeContainer::TemplateEdgeContainer(TemplateEdgeContainer const &other) : std::unordered_set<std::shared_ptr<TemplateEdge>>() {
             for (auto const& te : other) {
                 this->insert(std::make_shared<TemplateEdge>(*te));
             }
         }
+        
+        TemplateEdgeContainer& TemplateEdgeContainer::operator=(const TemplateEdgeContainer& other) {
+            this->clear();
+            for (auto const& te : other) {
+                this->insert(std::make_shared<TemplateEdge>(*te));
+            }
+            return *this;
+        }
     }
-}
+}

src/storm/storage/jani/TemplateEdgeContainer.h

@@ -11,6 +11,7 @@ namespace storm {
         struct TemplateEdgeContainer : public std::unordered_set<std::shared_ptr<TemplateEdge>>   {
             TemplateEdgeContainer() = default;
             TemplateEdgeContainer(TemplateEdgeContainer const& other);
+            TemplateEdgeContainer& operator=(TemplateEdgeContainer const& other);
         };
     }
 }

src/storm/utility/constants.cpp

@@ -55,12 +55,14 @@ namespace storm {
             return std::isnan(value);
         }
         
-        bool isAlmostZero(double const& a) {
-            return a < 1e-12 && a > -1e-12;
+        template<typename ValueType>
+        bool isAlmostZero(ValueType const& a) {
+            return a < convertNumber<ValueType>(1e-12) && a > -convertNumber<ValueType>(1e-12);
         }
         
-        bool isAlmostOne(double const& a) {
-            return a < (1.0 + 1e-12) && a > (1.0 - 1e-12);
+        template<typename ValueType>
+        bool isAlmostOne(ValueType const& a) {
+            return a < convertNumber<ValueType>(1.0 + 1e-12) && a > convertNumber<ValueType>(1.0 - 1e-12);
         }
         
         template<typename ValueType>
@@ -818,6 +820,16 @@ namespace storm {
             return std::move(value);
         }
         
+        template<>
+        bool isAlmostZero(storm::RationalFunction const& a) {
+            return a.isConstant() && isAlmostZero(convertNumber<RationalFunctionCoefficient>(a));
+        }
+        
+        template<>
+        bool isAlmostOne(storm::RationalFunction const& a) {
+            return a.isConstant() && isAlmostOne(convertNumber<RationalFunctionCoefficient>(a));
+        }
+        
         template<>
         std::pair<storm::RationalFunction, storm::RationalFunction> minmax(std::vector<storm::RationalFunction> const&) {
             STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Minimum/maximum for rational functions is not defined.");
@@ -885,6 +897,8 @@ namespace storm {
         template double infinity();
         template bool isOne(double const& value);
         template bool isZero(double const& value);
+        template bool isAlmostZero(double const& value);
+        template bool isAlmostOne(double const& value);
         template bool isConstant(double const& value);
         template bool isInfinity(double const& value);
         template bool isInteger(double const& number);
@@ -982,6 +996,8 @@ namespace storm {
         template bool isConstant(storm::ClnRationalNumber const& value);
         template bool isInfinity(storm::ClnRationalNumber const& value);
         template bool isNan(storm::ClnRationalNumber const& value);
+        template bool isAlmostZero(storm::ClnRationalNumber const& value);
+        template bool isAlmostOne(storm::ClnRationalNumber const& value);
         template storm::NumberTraits<ClnRationalNumber>::IntegerType convertNumber(storm::NumberTraits<ClnRationalNumber>::IntegerType const& number);
         template storm::ClnRationalNumber convertNumber(storm::ClnRationalNumber const& number);
         template storm::ClnRationalNumber simplify(storm::ClnRationalNumber value);
@@ -1008,6 +1024,8 @@ namespace storm {
         template bool isConstant(storm::GmpRationalNumber const& value);
         template bool isInfinity(storm::GmpRationalNumber const& value);
         template bool isNan(storm::GmpRationalNumber const& value);
+        template bool isAlmostZero(storm::GmpRationalNumber const& value);
+        template bool isAlmostOne(storm::GmpRationalNumber const& value);
         template storm::NumberTraits<GmpRationalNumber>::IntegerType convertNumber(storm::NumberTraits<GmpRationalNumber>::IntegerType const& number);
         template storm::GmpRationalNumber convertNumber(storm::GmpRationalNumber const& number);
         template storm::GmpRationalNumber simplify(storm::GmpRationalNumber value);

src/storm/utility/constants.h

@@ -88,9 +88,11 @@ namespace storm {
         template<typename ValueType>
         bool isNan(ValueType const& a);
 
-        bool isAlmostZero(double const& a);
+        template<typename ValueType>
+        bool isAlmostZero(ValueType const& a);
 
-        bool isAlmostOne(double const& a);
+        template<typename ValueType>
+        bool isAlmostOne(ValueType const& a);
         
         template<typename ValueType>
         bool isConstant(ValueType const& a);

src/storm/utility/shortestPaths.h

@@ -47,7 +47,7 @@ namespace storm {
                     if (predecessorNode != rhs.predecessorNode) {
                         return predecessorNode < rhs.predecessorNode;
                     }
-                    return predecessorK < predecessorK;
+                    return predecessorK < rhs.predecessorK;
                 }
 
                 bool operator==(const Path<T>& rhs) const {

version.cmake

@@ -1,4 +1,4 @@
 set(STORM_VERSION_MAJOR 1)
 set(STORM_VERSION_MINOR 6)
-set(STORM_VERSION_PATCH 0)
+set(STORM_VERSION_PATCH 2)