made sound value iteration work and respect the correct precision

src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp

@@ -13,9 +13,13 @@
 #include "storm/storage/memorystructure/MemoryStructureBuilder.h"
 #include "storm/storage/memorystructure/SparseModelMemoryProduct.h"
 #include "storm/storage/expressions/ExpressionManager.h"
+#include "storm/transformer/EndComponentEliminator.h"
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
 #include "storm/utility/graph.h"
+#include "storm/settings/SettingsManager.h"
+#include "storm/settings/modules/GeneralSettings.h"
+
 
 #include "storm/exceptions/InvalidPropertyException.h"
 #include "storm/exceptions/UnexpectedException.h"
@@ -383,10 +387,11 @@ namespace storm {
                 setReward0States(result, backwardTransitions);
                 checkRewardFiniteness(result, data.finiteRewardCheckObjectives, backwardTransitions);
                 
-                // Compute upper result bounds wherever this is necessarry
-                for (auto const& objIndex : data.upperResultBoundObjectives) {
-                    STORM_LOG_WARN("Computing upper result bound for an objective. TODO: this might not be necessary depending on the solution method.");
-                    result.objectives[objIndex].upperResultBound = computeUpperResultBound(result, objIndex, backwardTransitions);
+                // We compute upper result bounds if the 'sound' option has been enabled
+                if (storm::settings::getModule<storm::settings::modules::GeneralSettings>().isSoundSet()) {
+                    for (auto const& objIndex : data.upperResultBoundObjectives) {
+                        result.objectives[objIndex].upperResultBound = computeUpperResultBound(result, objIndex, backwardTransitions);
+                    }
                 }
                 
                 return result;
@@ -525,28 +530,40 @@ namespace storm {
                     auto nonZeroRewardStates = rewModel.getStatesWithZeroReward(transitions);
                     nonZeroRewardStates.complement();
                     auto expRewGreater0EStates = storm::utility::graph::performProbGreater0E(backwardTransitions, allStates, nonZeroRewardStates);
-                    auto expRew0AStates = ~expRewGreater0EStates;
+                    
+                    auto zeroRewardChoices = rewModel.getChoicesWithZeroReward(transitions);
+                    
+                    auto ecElimRes = storm::transformer::EndComponentEliminator<ValueType>::transform(transitions, expRewGreater0EStates, zeroRewardChoices, ~allStates);
+                    
+                    allStates.resize(ecElimRes.matrix.getRowGroupCount());
+                    storm::storage::BitVector outStates(allStates.size(), false);
+                    std::vector<ValueType> rew0StateProbs;
+                    rew0StateProbs.reserve(ecElimRes.matrix.getRowCount());
+                    for (uint64_t state = 0; state < allStates.size(); ++ state) {
+                        for (uint64_t choice = ecElimRes.matrix.getRowGroupIndices()[state]; choice < ecElimRes.matrix.getRowGroupIndices()[state + 1]; ++choice) {
+                            ValueType outProb = storm::utility::one<ValueType>() - ecElimRes.matrix.getRowSum(choice);
+                            if (!storm::utility::isZero(outProb)) {
+                                outStates.set(state, true);
+                            }
+                            rew0StateProbs.push_back(outProb);
+                        }
+                    }
                     
                     // An upper reward bound can only be computed if it is below infinity
-                    if (storm::utility::graph::performProb1A(transitions, transitions.getRowGroupIndices(), backwardTransitions, allStates, expRew0AStates).full()) {
-                        storm::storage::SparseMatrix<ValueType> submatrix = transitions.getSubmatrix(true, expRewGreater0EStates, expRewGreater0EStates);
-                        std::vector<ValueType> rew0StateProbs;
-                        rew0StateProbs.reserve(submatrix.getRowCount());
+                    if (storm::utility::graph::performProb1A(ecElimRes.matrix, ecElimRes.matrix.getRowGroupIndices(), ecElimRes.matrix.transpose(true), allStates, outStates).full()) {
+                        
+                        auto actionRewards = rewModel.getTotalRewardVector(transitions);
                         std::vector<ValueType> rewards;
-                        rewards.reserve(submatrix.getRowCount());
-                        for (auto const& state : expRewGreater0EStates) {
-                            for (uint64_t choice = transitions.getRowGroupIndices()[state]; choice < transitions.getRowGroupIndices()[state + 1]; ++choice) {
-                                rew0StateProbs.push_back(transitions.getConstrainedRowSum(choice, expRew0AStates));
-                                rewards.push_back(rewModel.getTotalStateActionReward(state, choice, transitions));
-                            }
+                        rewards.reserve(ecElimRes.matrix.getRowCount());
+                        for (auto row : ecElimRes.newToOldRowMapping) {
+                            rewards.push_back(actionRewards[row]);
                         }
                         
-                        storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ValueType> baier(submatrix, rewards, rew0StateProbs);
+                        storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ValueType> baier(ecElimRes.matrix, rewards, rew0StateProbs);
                         return baier.computeUpperBound();
                     }
                 }
                 
-                
                 // reaching this point means that we were not able to compute an upper bound
                 return boost::none;
 

src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp

@@ -12,6 +12,7 @@
 #include "storm/settings/SettingsManager.h"
 #include "storm/utility/export.h"
 #include "storm/settings/modules/IOSettings.h"
+#include "storm/settings/modules/GeneralSettings.h"
 
 #include "storm/exceptions/InvalidPropertyException.h"
 #include "storm/exceptions/InvalidOperationException.h"
@@ -48,8 +49,17 @@ namespace storm {
                 auto initEpoch = rewardUnfolding.getStartEpoch();
                 auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
                 EpochCheckingData cachedData;
+                ValueType precision = storm::utility::convertNumber<ValueType>(storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
+                if (storm::settings::getModule<storm::settings::modules::GeneralSettings>().isSoundSet()) {
+                    uint64_t denom = 0;
+                    for (uint64_t dim = 0; dim < rewardUnfolding.getEpochManager().getDimensionCount(); ++dim) {
+                        denom += rewardUnfolding.getEpochManager().getDimensionOfEpoch(initEpoch, dim) + 1;
+                    }
+                    precision = precision / storm::utility::convertNumber<ValueType>(denom);
+                }
+
                 for (auto const& epoch : epochOrder) {
-                    computeEpochSolution(epoch, weightVector, cachedData);
+                    computeEpochSolution(epoch, weightVector, cachedData, precision);
                     if (storm::settings::getModule<storm::settings::modules::IOSettings>().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) {
                         std::vector<ValueType> cdfEntry;
                         for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) {
@@ -84,7 +94,7 @@ namespace storm {
             }
             
             template <class SparseMdpModelType>
-            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData) {
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision) {
                 
                 ++numCheckedEpochs;
                 swEpochModelBuild.start();
@@ -169,7 +179,7 @@ namespace storm {
                 } else {
                     
                     
-                    updateCachedData(epochModel, cachedData, weightVector);
+                    updateCachedData(epochModel, cachedData, weightVector, precision);
                     
                     
                     // Formulate a min-max equation system max(A*x+b)=x for the weighted sum of the objectives
@@ -213,6 +223,7 @@ namespace storm {
                             subMatrix.convertToEquationSystem();
                         }
                         cachedData.linEqSolver = linEqSolverFactory.create(std::move(subMatrix));
+                        cachedData.linEqSolver->setPrecision(precision);
                         cachedData.linEqSolver->setCachingEnabled(true);
                     }
                     
@@ -277,7 +288,7 @@ namespace storm {
             }
 
             template <class SparseMdpModelType>
-            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector) {
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision) {
                 if (epochModel.epochMatrixChanged) {
                 
                     // Update the cached MinMaxSolver data
@@ -285,6 +296,7 @@ namespace storm {
                     cachedData.xMinMax.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
                     storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> minMaxSolverFactory;
                     cachedData.minMaxSolver = minMaxSolverFactory.create(epochModel.epochMatrix);
+                    cachedData.minMaxSolver->setPrecision(precision);
                     cachedData.minMaxSolver->setTrackScheduler(true);
                     cachedData.minMaxSolver->setCachingEnabled(true);
                     auto req = cachedData.minMaxSolver->getRequirements(storm::solver::EquationSystemType::StochasticShortestPath);

src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h

@@ -73,9 +73,9 @@ namespace storm {
                     std::vector<typename MultiDimensionalRewardUnfolding<ValueType, false>::SolutionType> solutions;
                 };
                 
-                void computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData);
+                void computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision);
                 
-                void updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector);
+                void updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision);
                 
                 storm::utility::Stopwatch swAll, swEpochModelBuild, swEpochModelAnalysis;
                 uint64_t numSchedChanges, numCheckedEpochs, numChecks;

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

@@ -24,6 +24,7 @@
  
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
+#include "storm/settings/modules/GeneralSettings.h"
 
 #include "storm/utility/Stopwatch.h"
 
@@ -87,6 +88,16 @@ namespace storm {
                 // initialize data that will be needed for each epoch
                 std::vector<ValueType> x, b;
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+
+                ValueType precision = storm::utility::convertNumber<ValueType>(storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
+                if (storm::settings::getModule<storm::settings::modules::GeneralSettings>().isSoundSet()) {
+                    uint64_t denom = 0;
+                    for (uint64_t dim = 0; dim < rewardUnfolding.getEpochManager().getDimensionCount(); ++dim) {
+                        denom += rewardUnfolding.getEpochManager().getDimensionOfEpoch(initEpoch, dim) + 1;
+                    }
+                    precision = precision / storm::utility::convertNumber<ValueType>(denom);
+                }
+
                 for (auto const& epoch : epochOrder) {
                     swBuild.start();
                     auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
@@ -132,6 +143,7 @@ namespace storm {
                         if (epochModel.epochMatrixChanged) {
                             x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
                             minMaxSolver = minMaxLinearEquationSolverFactory.create(epochModel.epochMatrix);
+                            minMaxSolver->setPrecision(precision);
                             minMaxSolver->setOptimizationDirection(dir);
                             minMaxSolver->setCachingEnabled(true);
                             minMaxSolver->setTrackScheduler(true);

src/storm/solver/AbstractEquationSolver.cpp

@@ -229,6 +229,13 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setPrecision(ValueType const& precision) {
+            STORM_LOG_DEBUG("Setting solver precision for a solver that does not support precisions.");
+        }
+
+        
+        
         template class AbstractEquationSolver<double>;
         template class AbstractEquationSolver<float>;
 

src/storm/solver/AbstractEquationSolver.h

@@ -155,6 +155,8 @@ namespace storm {
              * Shows progress if this solver is asked to do so.
              */
             void showProgressIterative(uint64_t iterations, boost::optional<uint64_t> const& bound = boost::none) const;
+            
+            virtual void setPrecision(ValueType const& precision);
 
         protected:
             /*!

src/storm/solver/GmmxxLinearEquationSolver.cpp

@@ -246,6 +246,11 @@ namespace storm {
             settings = newSettings;
         }
 
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::setPrecision(ValueType const& precision) {
+            settings.setPrecision(precision);
+        }
+        
         template<typename ValueType>
         GmmxxLinearEquationSolverSettings<ValueType> const& GmmxxLinearEquationSolver<ValueType>::getSettings() const {
             return settings;

src/storm/solver/GmmxxLinearEquationSolver.h

@@ -99,6 +99,8 @@ namespace storm {
             virtual void multiplyAndReduceGaussSeidel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) const override;
 
             void setSettings(GmmxxLinearEquationSolverSettings<ValueType> const& newSettings);
+            virtual void setPrecision(ValueType const& precision) override;
+
             GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;
 
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -245,6 +245,11 @@ namespace storm {
             return this->getSettings().getPrecision();
         }
 
+        template<typename ValueType>
+        void IterativeMinMaxLinearEquationSolver<ValueType>::setPrecision(ValueType const& precision) {
+            settings.setPrecision(precision);
+        }
+
         template<typename ValueType>
         bool IterativeMinMaxLinearEquationSolver<ValueType>::getRelative() const {
             return this->getSettings().getRelativeTerminationCriterion();

src/storm/solver/IterativeMinMaxLinearEquationSolver.h

@@ -56,6 +56,8 @@ namespace storm {
             virtual void clearCache() const override;
 
             ValueType getPrecision() const;
+            virtual void setPrecision(ValueType const& precision) override;
+
             bool getRelative() const;
             
             virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const override;

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -711,6 +711,11 @@ namespace storm {
         void NativeLinearEquationSolver<ValueType>::setSettings(NativeLinearEquationSolverSettings<ValueType> const& newSettings) {
             settings = newSettings;
         }
+
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::setPrecision(ValueType const& precision) {
+            settings.setPrecision(precision);
+        }
         
         template<typename ValueType>
         NativeLinearEquationSolverSettings<ValueType> const& NativeLinearEquationSolver<ValueType>::getSettings() const {

src/storm/solver/NativeLinearEquationSolver.h

@@ -66,6 +66,7 @@ namespace storm {
             
             void setSettings(NativeLinearEquationSolverSettings<ValueType> const& newSettings);
             NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
+            virtual void setPrecision(ValueType const& precision) override;
 
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
             virtual LinearEquationSolverRequirements getRequirements() const override;