added weight vector checker for reward bounded objectives

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

@@ -2,6 +2,7 @@
 
 #include "storm/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h"
 #include "storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h"
+#include "storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h"
 
 #include "storm/utility/macros.h"
 #include "storm/exceptions/NotSupportedException.h"
@@ -33,7 +34,11 @@ namespace storm {
             template <typename ModelType>
             template<typename VT, typename std::enable_if<std::is_same<ModelType, storm::models::sparse::Mdp<VT>>::value, int>::type>
             std::unique_ptr<PcaaWeightVectorChecker<ModelType>>  WeightVectorCheckerFactory<ModelType>::create(SparseMultiObjectivePreprocessorResult<ModelType> const& preprocessorResult) {
-                return std::make_unique<SparseMdpPcaaWeightVectorChecker<ModelType>>(preprocessorResult);
+                if (preprocessorResult.containsOnlyRewardObjectives()) {
+                    return std::make_unique<SparseMdpPcaaWeightVectorChecker<ModelType>>(preprocessorResult);
+                } else {
+                    return std::make_unique<SparseMdpRewardBoundedPcaaWeightVectorChecker<ModelType>>(preprocessorResult);
+                }
             }
             
             template <typename ModelType>

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

@@ -0,0 +1,129 @@
+#include "storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h"
+
+#include "storm/adapters/RationalFunctionAdapter.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/utility/macros.h"
+#include "storm/utility/vector.h"
+#include "storm/logic/Formulas.h"
+#include "storm/solver/MinMaxLinearEquationSolver.h"
+#include "storm/solver/LinearEquationSolver.h"
+
+
+#include "storm/exceptions/InvalidPropertyException.h"
+#include "storm/exceptions/InvalidOperationException.h"
+#include "storm/exceptions/IllegalArgumentException.h"
+#include "storm/exceptions/NotSupportedException.h"
+#include "storm/exceptions/UnexpectedException.h"
+
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template <class SparseMdpModelType>
+            SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::SparseMdpRewardBoundedPcaaWeightVectorChecker(SparseMultiObjectivePreprocessorResult<SparseMdpModelType> const& preprocessorResult) : PcaaWeightVectorChecker<SparseMdpModelType>(preprocessorResult.objectives), rewardUnfolding(*preprocessorResult.preprocessedModel, this->objectives, storm::storage::BitVector(preprocessorResult.preprocessedModel->getNumberOfChoices(), true), preprocessorResult.reward0EStates) {
+            
+            }
+            
+            template <class SparseMdpModelType>
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::check(std::vector<ValueType> const& weightVector) {
+                auto initEpoch = rewardUnfolding.getStartEpoch();
+                auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
+                for (auto const& epoch : epochOrder) {
+                    computeEpochSolution(epoch, weightVector);
+                }
+                
+                auto solution = rewardUnfolding.getInitialStateResult(initEpoch);
+                // Todo: we currently assume precise results...
+                underApproxResult = solution.objectiveValues;
+                overApproxResult = solution.objectiveValues;
+                
+            }
+            
+            template <class SparseMdpModelType>
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType>::Epoch const& epoch, std::vector<ValueType> const& weightVector) {
+                auto const& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
+                std::vector<typename MultiDimensionalRewardUnfolding<ValueType>::SolutionType> result(epochModel.epochMatrix.getRowGroupCount());
+                
+                
+                // Formulate a min-max equation system max(A*x+b)=x for the weighted sum of the objectives
+                std::vector<ValueType> b(epochModel.epochMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                    ValueType weight = storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType()) ? -weightVector[objIndex] : weightVector[objIndex];
+                    if (!storm::utility::isZero(weight)) {
+                        std::vector<ValueType> const& objectiveReward = epochModel.objectiveRewards[objIndex];
+                        for (auto const& choice : epochModel.objectiveRewardFilter[objIndex]) {
+                            b[choice] += weight * objectiveReward[choice];
+                        }
+                    }
+                }
+                auto stepSolutionIt = epochModel.stepSolutions.begin();
+                for (auto const& choice : epochModel.stepChoices) {
+                    b[choice] += stepSolutionIt->weightedValue;
+                    ++stepSolutionIt;
+                }
+                
+                // Invoke the min max solver
+                storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> minMaxSolverFactory;
+                auto minMaxSolver = minMaxSolverFactory.create(epochModel.epochMatrix);
+                minMaxSolver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize);
+                minMaxSolver->setTrackScheduler(true);
+                //minMaxSolver->setCachingEnabled(true);
+                std::vector<ValueType> x(result.size(), storm::utility::zero<ValueType>());
+                minMaxSolver->solveEquations(x, b);
+                for (uint64_t state = 0; state < x.size(); ++state) {
+                    result[state].weightedValue = x[state];
+                }
+                
+                // Formulate for each objective the linear equation system induced by the performed choices
+                auto const& choices = minMaxSolver->getSchedulerChoices();
+                storm::storage::SparseMatrix<ValueType> subMatrix = epochModel.epochMatrix.selectRowsFromRowGroups(choices, true);
+                subMatrix.convertToEquationSystem();
+                storm::solver::GeneralLinearEquationSolverFactory<ValueType> linEqSolverFactory;
+                auto linEqSolver = linEqSolverFactory.create(std::move(subMatrix));
+                b.resize(choices.size());
+                // TODO: start with a better initial guess
+                x.resize(choices.size());
+                for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                    std::vector<ValueType> const& objectiveReward = epochModel.objectiveRewards[objIndex];
+                    for (uint64_t state = 0; state < choices.size(); ++state) {
+                        uint64_t choice = epochModel.epochMatrix.getRowGroupIndices()[state] + choices[state];
+                        if (epochModel.objectiveRewardFilter[objIndex].get(choice)) {
+                            b[state] = objectiveReward[choice];
+                        } else {
+                            b[state] = storm::utility::zero<ValueType>();
+                        }
+                        if (epochModel.stepChoices.get(choice)) {
+                            b[state] += epochModel.stepSolutions[epochModel.stepChoices.getNumberOfSetBitsBeforeIndex(choice)].objectiveValues[objIndex];
+                        }
+                    }
+                    linEqSolver->solveEquations(x, b);
+                    for (uint64_t state = 0; state < choices.size(); ++state) {
+                        result[state].objectiveValues.push_back(x[state]);
+                    }
+                }
+                
+                rewardUnfolding.setSolutionForCurrentEpoch(result);
+            }
+
+            template <class SparseMdpModelType>
+            std::vector<typename SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::ValueType> SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::getUnderApproximationOfInitialStateResults() const {
+                STORM_LOG_THROW(underApproxResult, storm::exceptions::InvalidOperationException, "Tried to retrieve results but check(..) has not been called before.");
+                return underApproxResult.get();
+            }
+            
+            template <class SparseMdpModelType>
+            std::vector<typename SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::ValueType> SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::getOverApproximationOfInitialStateResults() const {
+                STORM_LOG_THROW(overApproxResult, storm::exceptions::InvalidOperationException, "Tried to retrieve results but check(..) has not been called before.");
+                return overApproxResult.get();
+            }
+            
+            template class SparseMdpRewardBoundedPcaaWeightVectorChecker<storm::models::sparse::Mdp<double>>;
+#ifdef STORM_HAVE_CARL
+            template class SparseMdpRewardBoundedPcaaWeightVectorChecker<storm::models::sparse::Mdp<storm::RationalNumber>>;
+#endif
+        
+        }
+    }
+}

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

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <vector>
+
+#include "storm/modelchecker/multiobjective/pcaa/PcaaWeightVectorChecker.h"
+#include "storm/modelchecker/multiobjective/rewardbounded/MultiDimensionalRewardUnfolding.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            /*!
+             * Helper Class that takes preprocessed Pcaa data and a weight vector and ...
+             * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives
+             * - extracts the scheduler that induces this maximum
+             * - computes for each objective the value induced by this scheduler
+             */
+            template <class SparseMdpModelType>
+            class SparseMdpRewardBoundedPcaaWeightVectorChecker : public PcaaWeightVectorChecker<SparseMdpModelType> {
+            public:
+                typedef typename SparseMdpModelType::ValueType ValueType;
+                typedef typename SparseMdpModelType::RewardModelType RewardModelType;
+            
+                SparseMdpRewardBoundedPcaaWeightVectorChecker(SparseMultiObjectivePreprocessorResult<SparseMdpModelType> const& preprocessorResult);
+
+                virtual ~SparseMdpRewardBoundedPcaaWeightVectorChecker() = default;
+
+                /*!
+                 * - computes the optimal expected reward w.r.t. the weighted sum of the rewards of the individual objectives
+                 * - extracts the scheduler that induces this optimum
+                 * - computes for each objective the value induced by this scheduler
+                 */
+                virtual void check(std::vector<ValueType> const& weightVector) override;
+                
+                /*!
+                 * Retrieves the results of the individual objectives at the initial state of the given model.
+                 * Note that check(..) has to be called before retrieving results. Otherwise, an exception is thrown.
+                 * Also note that there is no guarantee that the under/over approximation is in fact correct
+                 * as long as the underlying solution methods are unsound (e.g., standard value iteration).
+                 */
+                virtual std::vector<ValueType> getUnderApproximationOfInitialStateResults() const override;
+                virtual std::vector<ValueType> getOverApproximationOfInitialStateResults() const override;
+                
+            private:
+                
+                void computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType>::Epoch const& epoch, std::vector<ValueType> const& weightVector);
+                
+                
+                MultiDimensionalRewardUnfolding<ValueType> rewardUnfolding;
+                
+                boost::optional<std::vector<ValueType>> underApproxResult;
+                boost::optional<std::vector<ValueType>> overApproxResult;
+
+                
+            };
+            
+        }
+    }
+}