better lower/upper result bounds

src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp

@@ -530,46 +530,6 @@ namespace storm {
                     auto const& rewModel = result.preprocessedModel->getRewardModel(result.objectives[objIndex].formula->asRewardOperatorFormula().getRewardModelName());
                     auto actionRewards = rewModel.getTotalRewardVector(transitions);
                     
-                    // TODO: Consider cumulative reward formulas less naively
-                    if (result.objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) {
-                        auto const& cumulativeRewardFormula = result.objectives[objIndex].formula->getSubformula().asCumulativeRewardFormula();
-                        ValueType rewardBound = cumulativeRewardFormula.template getBound<ValueType>();
-                        if (cumulativeRewardFormula.getTimeBoundReference().isRewardBound()) {
-                            auto const& costModel = result.preprocessedModel->getRewardModel(cumulativeRewardFormula.getTimeBoundReference().getRewardName());
-                            if (!costModel.hasTransitionRewards()) {
-                                auto actionCosts = costModel.getTotalRewardVector(transitions);
-                                typename  SparseModelType::ValueType largestRewardPerCost = storm::utility::zero<typename SparseModelType::ValueType>();
-                                bool isFinite = true;
-                                for (auto rewIt = actionRewards.begin(), costIt = actionCosts.begin(); rewIt != actionRewards.end(); ++rewIt, ++costIt) {
-                                    if (!storm::utility::isZero(*rewIt)) {
-                                        if (storm::utility::isZero(*costIt)) {
-                                            isFinite = false;
-                                            break;
-                                        }
-                                        ValueType rewardPerCost = *rewIt / *costIt;
-                                        largestRewardPerCost = std::max(largestRewardPerCost, rewardPerCost);
-                                    }
-                                }
-                                
-                                if (isFinite) {
-                                    ValueType newResultBound = largestRewardPerCost * rewardBound;
-                                    if (upperBound) {
-                                        upperBound = std::min(upperBound.get(), newResultBound);
-                                    } else {
-                                        upperBound = newResultBound;
-                                    }
-                                }
-                            }
-                        } else {
-                            ValueType newResultBound = (*std::max_element(actionRewards.begin(), actionRewards.end())) * rewardBound;
-                            if (upperBound) {
-                                upperBound = std::min(upperBound.get(), newResultBound);
-                            } else {
-                                upperBound = newResultBound;
-                            }
-                        }
-                    }
-                    
                     if (result.objectives[objIndex].formula->getSubformula().isTotalRewardFormula() || result.objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) {
                         // We have to eliminate ECs here to treat zero-reward ECs
  
@@ -590,11 +550,21 @@ namespace storm {
                         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);
+                                // Check whether the choice lead to a state with expRew 0 in the original model
+                                bool isOutChoice = false;
+                                uint64_t originalModelChoice = ecElimRes.newToOldRowMapping[choice];
+                                for (auto const& entry : transitions.getRow(originalModelChoice)) {
+                                    if (!expRewGreater0EStates.get(entry.getColumn())) {
+                                        isOutChoice = true;
+                                        outStates.set(state, true);
+                                        rew0StateProbs.push_back(storm::utility::one<ValueType>() - ecElimRes.matrix.getRowSum(choice));
+                                        assert (!storm::utility::isZero(rew0StateProbs.back()));
+                                        break;
+                                    }
+                                }
+                                if (!isOutChoice) {
+                                    rew0StateProbs.push_back(storm::utility::zero<ValueType>());
                                 }
-                                rew0StateProbs.push_back(outProb);
                             }
                         }
                         

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

@@ -34,6 +34,12 @@ namespace storm {
                 STORM_LOG_THROW(preprocessorResult.rewardFinitenessType == SparseMultiObjectivePreprocessorResult<SparseMdpModelType>::RewardFinitenessType::AllFinite, storm::exceptions::NotSupportedException, "There is a scheduler that yields infinite reward for one  objective. This is not supported.");
                 STORM_LOG_THROW(preprocessorResult.preprocessedModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "The model has multiple initial states.");
                 
+                // Update the objective bounds with what the reward unfolding can compute
+                for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                    this->objectives[objIndex].lowerResultBound = rewardUnfolding.getLowerObjectiveBound(objIndex);
+                    this->objectives[objIndex].upperResultBound = rewardUnfolding.getUpperObjectiveBound(objIndex);
+                }
+                
                 numCheckedEpochs = 0;
                 numChecks = 0;
             }

src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp

@@ -67,7 +67,7 @@ namespace storm {
                 }
                 auto formula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(checkTask.getFormula().asSharedPointer(), opInfo);
                 helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
-                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(helper::SolutionType::UntilProbabilities, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
+                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
                 return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
             } else {
                 STORM_LOG_THROW(!pathFormula.hasLowerBound() && pathFormula.hasUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have single upper time bound.");
@@ -146,7 +146,7 @@ namespace storm {
                 }
                 auto formula = std::make_shared<storm::logic::RewardOperatorFormula>(checkTask.getFormula().asSharedPointer(), checkTask.getRewardModel(), opInfo);
                 helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
-                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(helper::SolutionType::ExpectedRewards, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
+                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
                 return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
             } else {
                 

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

@@ -134,7 +134,7 @@ namespace storm {
             }
             
             template<typename ValueType>
-            std::vector<ValueType> analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ValueType const& precision, SolutionType const& type) {
+            std::vector<ValueType> analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ValueType const& precision, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
  
                 // Update some data for the case that the Matrix has changed
                 if (epochModel.epochMatrixChanged) {
@@ -146,14 +146,13 @@ namespace storm {
                     minMaxSolver->setCachingEnabled(true);
                     minMaxSolver->setTrackScheduler(true);
                     auto req = minMaxSolver->getRequirements(dir);
-                    minMaxSolver->setLowerBound(storm::utility::zero<ValueType>());
-                    req.clearLowerBounds();
-                    if (type == SolutionType::UntilProbabilities) {
-                        minMaxSolver->setUpperBound(storm::utility::one<ValueType>());
+                    if (lowerBound) {
+                        minMaxSolver->setLowerBound(lowerBound.get());
+                        req.clearLowerBounds();
+                    }
+                    if (upperBound) {
+                        minMaxSolver->setUpperBound(upperBound.get());
                         req.clearUpperBounds();
-                    } else if (type == SolutionType::ExpectedRewards) {
-                        // TODO
-                        STORM_LOG_WARN_COND(!req.requiresUpperBounds(), "Upper bounds for expected reward are not specified.");
                     }
                     STORM_LOG_THROW(req.empty(), storm::exceptions::UncheckedRequirementException, "At least one requirement was not checked.");
                     minMaxSolver->setRequirementsChecked();
@@ -182,8 +181,14 @@ namespace storm {
             }
             
             template<typename ValueType>
-            std::map<storm::storage::sparse::state_type, ValueType> SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+            std::map<storm::storage::sparse::state_type, ValueType> SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
                 storm::utility::Stopwatch swAll(true), swBuild, swCheck;
+                
+                // Get lower and upper bounds for the solution.
+                auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+                auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+                
+                // Initialize epoch models
                 auto initEpoch = rewardUnfolding.getStartEpoch();
                 auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
                 
@@ -205,7 +210,7 @@ namespace storm {
                     if (epochModel.epochMatrix.getEntryCount() == 0) {
                         rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialEpochModel<ValueType>(dir, epochModel));
                     } else {
-                        rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialEpochModel<ValueType>(dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, precision, type));
+                        rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialEpochModel<ValueType>(dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, precision, lowerBound, upperBound));
                     }
                     swCheck.stop();
                     ++numCheckedEpochs;

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

@@ -38,7 +38,7 @@ namespace storm {
                 
                 static std::vector<ValueType> computeStepBoundedUntilProbabilities(storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint());
 
-                static std::map<storm::storage::sparse::state_type, ValueType> computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                static std::map<storm::storage::sparse::state_type, ValueType> computeRewardBoundedValues(OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
                 
                 static std::vector<ValueType> computeNextProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
 

src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp

@@ -9,6 +9,8 @@
 
 #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h"
 #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h"
+
 #include "storm/transformer/EndComponentEliminator.h"
 
 #include "storm/exceptions/UnexpectedException.h"
@@ -548,6 +550,103 @@ namespace storm {
                     return precision / storm::utility::convertNumber<ValueType>(sumOfDimensions);
                 }
                 
+                template<typename ValueType, bool SingleObjectiveMode>
+                boost::optional<ValueType> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getUpperObjectiveBound(uint64_t objectiveIndex) {
+                    auto& objective = this->objectives[objectiveIndex];
+                    if (!objective.upperResultBound) {
+                        if (objective.formula->isProbabilityOperatorFormula()) {
+                            objective.upperResultBound = storm::utility::one<ValueType>();
+                        } else if (objective.formula->isRewardOperatorFormula()) {
+                            auto const& rewModel = this->model.getRewardModel(objective.formula->asRewardOperatorFormula().getRewardModelName());
+                            auto actionRewards = rewModel.getTotalRewardVector(this->model.getTransitionMatrix());
+                            if (objective.formula->getSubformula().isCumulativeRewardFormula()) {
+                                // Try to get an upper bound by computing the maximal reward achievable within one epoch step
+                                auto const& cumulativeRewardFormula = objective.formula->getSubformula().asCumulativeRewardFormula();
+                                ValueType rewardBound = cumulativeRewardFormula.template getBound<ValueType>();
+                                if (cumulativeRewardFormula.getTimeBoundReference().isRewardBound()) {
+                                    auto const& costModel = this->model.getRewardModel(cumulativeRewardFormula.getTimeBoundReference().getRewardName());
+                                    if (!costModel.hasTransitionRewards()) {
+                                        auto actionCosts = costModel.getTotalRewardVector(this->model.getTransitionMatrix());
+                                        ValueType largestRewardPerCost = storm::utility::zero<ValueType>();
+                                        bool isFinite = true;
+                                        for (auto rewIt = actionRewards.begin(), costIt = actionCosts.begin(); rewIt != actionRewards.end(); ++rewIt, ++costIt) {
+                                            if (!storm::utility::isZero(*rewIt)) {
+                                                if (storm::utility::isZero(*costIt)) {
+                                                    isFinite = false;
+                                                    break;
+                                                }
+                                                ValueType rewardPerCost = *rewIt / *costIt;
+                                                largestRewardPerCost = std::max(largestRewardPerCost, rewardPerCost);
+                                            }
+                                        }
+                                        if (isFinite) {
+                                            objective.upperResultBound = largestRewardPerCost * rewardBound;
+                                        }
+                                    }
+                                } else {
+                                    objective.upperResultBound = (*std::max_element(actionRewards.begin(), actionRewards.end())) * rewardBound;
+                                }
+                                
+                                // If we could not find an upper bound, try to get an upper bound for the unbounded case
+                                if (!objective.upperResultBound) {
+                                    storm::storage::BitVector allStates(model.getNumberOfStates(), true);
+                                    // Get the set of states from which reward is reachable
+                                    auto nonZeroRewardStates = rewModel.getStatesWithZeroReward(model.getTransitionMatrix());
+                                    nonZeroRewardStates.complement();
+                                    auto expRewGreater0EStates = storm::utility::graph::performProbGreater0E(model.getBackwardTransitions(), allStates, nonZeroRewardStates);
+                                    // Eliminate zero-reward ECs
+                                    auto zeroRewardChoices = rewModel.getChoicesWithZeroReward(model.getTransitionMatrix());
+                                    auto ecElimRes = storm::transformer::EndComponentEliminator<ValueType>::transform(model.getTransitionMatrix(), 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) {
+                                            // Check whether the choice lead to a state with expRew 0 in the original model
+                                            bool isOutChoice = false;
+                                            uint64_t originalModelChoice = ecElimRes.newToOldRowMapping[choice];
+                                            for (auto const& entry : model.getTransitionMatrix().getRow(originalModelChoice)) {
+                                                if (!expRewGreater0EStates.get(entry.getColumn())) {
+                                                    isOutChoice = true;
+                                                    outStates.set(state, true);
+                                                    rew0StateProbs.push_back(storm::utility::one<ValueType>() - ecElimRes.matrix.getRowSum(choice));
+                                                    assert (!storm::utility::isZero(rew0StateProbs.back()));
+                                                    break;
+                                                }
+                                            }
+                                            if (!isOutChoice) {
+                                                rew0StateProbs.push_back(storm::utility::zero<ValueType>());
+                                            }
+                                        }
+                                    }
+                                    // An upper reward bound can only be computed if it is below infinity
+                                    if (storm::utility::graph::performProb1A(ecElimRes.matrix, ecElimRes.matrix.getRowGroupIndices(), ecElimRes.matrix.transpose(true), allStates, outStates).full()) {
+                                        std::vector<ValueType> rewards;
+                                        rewards.reserve(ecElimRes.matrix.getRowCount());
+                                        for (auto row : ecElimRes.newToOldRowMapping) {
+                                            rewards.push_back(actionRewards[row]);
+                                        }
+                                        storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ValueType> baier(ecElimRes.matrix, rewards, rew0StateProbs);
+                                        objective.upperResultBound = baier.computeUpperBound();
+                                    }
+                                }
+                            }
+                        }
+                    }
+                    return objective.upperResultBound;
+                }
+
+                template<typename ValueType, bool SingleObjectiveMode>
+                boost::optional<ValueType> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getLowerObjectiveBound(uint64_t objectiveIndex) {
+                    auto& objective = this->objectives[objectiveIndex];
+                    if (!objective.lowerResultBound) {
+                        objective.lowerResultBound = storm::utility::zero<ValueType>();
+                    }
+                    return objective.lowerResultBound;
+                }
+
+                
                 template<typename ValueType, bool SingleObjectiveMode>
                 void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions) {
                     STORM_LOG_ASSERT(currentEpoch, "Tried to set a solution for the current epoch, but no epoch was specified before.");

src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h

@@ -58,6 +58,12 @@ namespace storm {
                      */
                     ValueType getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision);
                     
+                    /*!
+                     * Returns an upper/lower bound for the objective result in every state (if this bound could be computed)
+                     */
+                    boost::optional<ValueType> getUpperObjectiveBound(uint64_t objectiveIndex = 0);
+                    boost::optional<ValueType> getLowerObjectiveBound(uint64_t objectiveIndex = 0);
+                    
                     void setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions);
                     SolutionType const& getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
                     SolutionType const& getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);