From 6aeb75e3bddfcbe73ced9bdbac47cdadfcedb45a Mon Sep 17 00:00:00 2001
From: TimQu <tim.quatmann@cs.rwth-aachen.de>
Date: Tue, 5 Feb 2019 19:14:15 +0100
Subject: [PATCH] quantiles: Supporting two-dimensional quantiles with the same
 optimization direction of quantile bounds (max,max or min,min).

---
 .../MultiDimensionalRewardUnfolding.cpp       |   4 +-
 .../helper/rewardbounded/ProductModel.cpp     |   2 +-
 .../helper/rewardbounded/QuantileHelper.cpp   | 288 +++++++++++++++---
 .../helper/rewardbounded/QuantileHelper.h     |   1 +
 4 files changed, 255 insertions(+), 40 deletions(-)

diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
index c710d1a04..1d57dd2b7 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
@@ -230,10 +230,10 @@ namespace storm {
                         }
                         
                         if (!bound.containsVariables()) {
-                            ValueType discretizedBound = storm::utility::convertNumber<ValueType>(bound.evaluateAsRational());
                             // We always consider upper bounds to be non-strict and lower bounds to be strict.
-                            // Thus, >=N would become >N-1. However, note that the case N=0 needs extra treatment
+                            // Thus, >=N would become >N-1. However, note that the case N=0 is treated separately.
                             if (dimensions[dim].isBounded) {
+                                ValueType discretizedBound = storm::utility::convertNumber<ValueType>(bound.evaluateAsRational());
                                 discretizedBound /= dimensions[dim].scalingFactor;
                                 if (storm::utility::isInteger(discretizedBound)) {
                                     if (isStrict == dimensions[dim].isUpperBounded) {
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
index a2c12698b..57fbf4980 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
@@ -511,7 +511,7 @@ namespace storm {
                     for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
                         if (epochManager.isBottomDimensionEpochClass(epochClass, dim)) {
                             bottomDimensions.set(dim, true);
-                            if (dimensions[dim].isBounded) {
+                            if (dimensions[dim].isBounded && dimensions[dim].maxValue) {
                                 considerInitialStates = false;
                             }
                         }
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index d20f4a7da..a06907e2f 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -4,7 +4,6 @@
 #include <vector>
 #include <memory>
 #include <boost/optional.hpp>
-#include <storm/exceptions/NotSupportedException.h>
 
 #include "storm/models/sparse/Mdp.h"
 #include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
@@ -17,6 +16,8 @@
 #include "storm/logic/ProbabilityOperatorFormula.h"
 #include "storm/logic/BoundedUntilFormula.h"
 
+#include "storm/exceptions/NotSupportedException.h"
+
 namespace storm {
     namespace modelchecker {
         namespace helper {
@@ -25,10 +26,20 @@ namespace storm {
                 template<typename ModelType>
                 QuantileHelper<ModelType>::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) {
                     // Intentionally left empty
-                    /* Todo: Current assumption:
+                    STORM_LOG_THROW(quantileFormula.getSubformula().isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs probability operator inside. The formula " << quantileFormula << " is not supported.");
+                    auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
+                    STORM_LOG_THROW(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, storm::exceptions::NotSupportedException, "Probability operator inside quantile formula needs to have bound > or <=.");
+                    STORM_LOG_THROW(probOpFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs bounded until probability operator formula as subformula. The formula " << quantileFormula << " is not supported.");
+                    auto const& boundedUntilFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
+                    
+                    // Only > and <= are supported for upper bounds. This is to make sure that Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
+                    // Only >= and < are supported for lower bounds. (EC construction..)
+                    // TODO
+                    
+                    /* Todo: Current assumptions:
                      * Subformula is always prob operator with bounded until
                      * Each bound variable occurs at most once (change this?)
-                     * cost bounds can either be <= or >  (the epochs returned by the reward unfolding assume this. One could translate them, though)
+                     * cost bounds are assumed to be always non-strict (the epochs returned by the reward unfolding assume strict lower and non-strict upper cost bounds, though...)
                      * 'reasonable' quantile (e.g. not quantile(max B, Pmax>0.5 [F <=B G]) (we could also filter out these things early on)
                      */
                 }
@@ -141,10 +152,8 @@ namespace storm {
                             result = {{storm::utility::convertNumber<ValueType>(quantileRes.first) * quantileRes.second}};
                         } else if (maxProbSatisfiesFormula) {
                             // i.e., all bound values satisfy the formula
-                            bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
-                            bool upperCostBound = boundedUntilOperator.getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
-                            if (minimizingRewardBound) {
-                                result = {{ (upperCostBound ? storm::utility::zero<ValueType>() : -storm::utility::one<ValueType>())}};
+                            if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) {
+                                result = {{ storm::utility::zero<ValueType>() }};
                             } else {
                                 result = {{ storm::utility::infinity<ValueType>()}};
                             }
@@ -153,52 +162,246 @@ namespace storm {
                             result = {{}};
                         }
                     } else if (getOpenDimensions().getNumberOfSetBits() == 2) {
-                    
+                        result = computeTwoDimensionalQuantile(env);
                     } else {
                         STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The quantile formula considers " << getOpenDimensions().getNumberOfSetBits() << " open dimensions. Only one- or two-dimensional quantiles are supported.");
                     }
                     return result;
-                    /*
-                                        return unboundedFormula->getBound().isSatisfied(numericResult);
-                    
-                    if (!checkUnboundedValue(env)) {
-                        STORM_LOG_INFO("Bound is not satisfiable.");
-                        result.emplace_back();
-                        for (auto const& bv : quantileFormula.getBoundVariables()) {
-                            result.front().push_back(storm::solver::minimize(bv.first) ? storm::utility::infinity<ValueType>() : -storm::utility::infinity<ValueType>());
+                }
+
+                template<typename ValueType>
+                void filterDominatedPoints(std::vector<std::vector<ValueType>>& points, std::vector<storm::solver::OptimizationDirection> const& dirs) {
+                    std::vector<std::vector<ValueType>> result;
+                    // Note: this is slow and not inplace but also most likely not performance critical
+                    for (auto const& p1 : points) {
+                        bool p1Dominated = false;
+                        for (auto const& p2 : points) {
+                            assert(p1.size() == p2.size());
+                            bool p2DominatesP1 = false;
+                            for (uint64_t i = 0; i < dirs.size(); ++i) {
+                                if (storm::solver::minimize(dirs[i]) ? p2[i] <= p1[i] : p2[i] >= p1[i]) {
+                                    if (p2[i] != p1[i]) {
+                                        p2DominatesP1 = true;
+                                    }
+                                } else {
+                                    p2DominatesP1 = false;
+                                    break;
+                                }
+                            }
+                            if (p2DominatesP1) {
+                                p1Dominated = true;
+                                break;
+                            }
                         }
-                    } else {
-                        std::vector<bool> limitProbCheckResults;
-                        for (uint64_t dim = 0; dim < getDimension(); ++dim) {
-                            storm::storage::BitVector dimAsBitVector(getDimension(), false);
-                            dimAsBitVector.set(dim, true);
-                            limitProbCheckResults.push_back(checkLimitProbability(env, dimAsBitVector));
-                            auto quantileRes = computeQuantileForDimension(env, dim);
-                            std::cout << "Quantile for dim " << dim << " is " << (storm::utility::convertNumber<ValueType>(quantileRes.first) * quantileRes.second) << std::endl;
+                        if (!p1Dominated) {
+                            result.push_back(p1);
                         }
-                        
-                        result = {{27}};
                     }
-                    return result;*/
+                    points = std::move(result);
                 }
+                
+                template<typename ModelType>
+                std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment const& env) {
+                    std::vector<std::vector<ValueType>> result;
 
+                    auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
+                    storm::logic::ComparisonType probabilityBound = probOpFormula.getBound().comparisonType;
+                    auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
+                    std::vector<storm::storage::BitVector> dimensionsAsBitVector;
+                    std::vector<uint64_t> dimensions;
+                    std::vector<storm::solver::OptimizationDirection> optimizationDirections;
+                    std::vector<bool> lowerCostBounds;
+                    for (auto const& dirVar : quantileFormula.getBoundVariables()) {
+                        dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
+                        STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
+                        dimensions.push_back(*dimensionsAsBitVector.back().begin());
+                        lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
+                        optimizationDirections.push_back(dirVar.first);
+                    }
+                    STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
+                    if (optimizationDirections[0] == optimizationDirections[1]) {
+                        if (lowerCostBounds[0] == lowerCostBounds[1]) {
+                            // TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0])
+                            
+                            bool maxmaxProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, storm::storage::BitVector(getDimension(), false)));
+                            bool minminProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1]));
+                            if (maxmaxProbSatisfiesFormula != minminProbSatisfiesFormula) {
+                                std::vector<std::pair<int64_t, typename ModelType::ValueType>> singleQuantileValues;
+                                std::vector<std::vector<int64_t>> resultPoints;
+                                const int64_t infinity = std::numeric_limits<int64_t>().max(); // use this value to represent infinity in a result point
+                                for (uint64_t i = 0; i < 2; ++i) {
+                                    // find out whether the bounds B_i = 0 and B_1-i = infinity satisfy the formula
+                                    uint64_t indexToMinimizeProb = lowerCostBounds[i] ? (1-i) : i;
+                                    bool zeroInfSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, dimensionsAsBitVector[indexToMinimizeProb]));
+                                    std::cout << "Formula sat is " << zeroInfSatisfiesFormula << " and lower bound is " << storm::logic::isLowerBound(probabilityBound) << std::endl;
+                                    if (zeroInfSatisfiesFormula == storm::solver::minimize(optimizationDirections[0])) {
+                                        // There is bound value b such that the point B_i=0 and B_1-i = b is part of the result
+                                        singleQuantileValues.emplace_back(0, storm::utility::zero<ValueType>());
+                                    } else {
+                                        // Compute quantile where 1-i approaches infinity
+                                        std::cout << "Computing quantile for single dimension " << dimensions[i] << std::endl;
+                                        singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]));
+                                        std::cout << ".. Result is " << singleQuantileValues.back().first << std::endl;
+                                        // When maximizing bounds, the computed quantile value is sat for all values of the other bound.
+                                        if (!storm::solver::minimize(optimizationDirections[i])) {
+                                            std::vector<int64_t> newResultPoint(2);
+                                            newResultPoint[i] = singleQuantileValues.back().first;
+                                            newResultPoint[1-i] = infinity;
+                                            resultPoints.push_back(newResultPoint);
+                                            // Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i]
+                                            ++singleQuantileValues.back().first;
+                                        }
+                                    }
+                                    // Decrease the value for lower cost bounds to convert >= to >
+                                    if (lowerCostBounds[i]) {
+                                        --singleQuantileValues[i].first;
+                                    }
+                                }
+                                std::cout << "Found starting point. Beginning 2D exploration" << std::endl;
+                                
+                                MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None)));
+                        
+                                // initialize data that will be needed for each epoch
+                                auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+                                auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+                                std::vector<ValueType> x, b;
+                                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+                    
+                                std::vector<int64_t> epochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0).
+                                std::cout << "Starting quantile exploration with: (" << epochValues[0] << ", " << epochValues[1] << ")." << std::endl;
+                                std::set<typename EpochManager::Epoch> exploredEpochs;
+                                while (true) {
+                                    auto currEpoch = rewardUnfolding.getStartEpoch(true);
+                                    for (uint64_t i = 0; i < 2; ++i) {
+                                        if (epochValues[i] >= 0) {
+                                            rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) epochValues[i]);
+                                        } else {
+                                            rewardUnfolding.getEpochManager().setBottomDimension(currEpoch, dimensions[i]);
+                                        }
+                                    }
+                                    auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
+                                    for (auto const& epoch : epochOrder) {
+                                        if (exploredEpochs.count(epoch) == 0) {
+                                            auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
+                                            rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+                                            exploredEpochs.insert(epoch);
+                                        }
+                                    }
+                                    ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
+                                    std::cout << "Numeric result for epoch " << rewardUnfolding.getEpochManager().toString(currEpoch) << " is " << currValue << std::endl;
+                                    bool propertySatisfied = probOpFormula.getBound().isSatisfied(currValue);
+                                    // If the reward bounds should be as small as possible, we should stop as soon as the property is satisfied.
+                                    // If the reward bounds should be as large as possible, we should stop as soon as the property is violated and then go a step backwards
+                                    if (storm::solver::minimize(optimizationDirections[0]) == propertySatisfied) {
+                                        // We found another point for the result!
+                                        resultPoints.push_back(epochValues);
+                                        std::cout << "Found another result point: (" << epochValues[0] << ", " << epochValues[1] << ")." << std::endl;
+                                        if (epochValues[1] == singleQuantileValues[1].first) {
+                                            break;
+                                        } else {
+                                            ++epochValues[0];
+                                            epochValues[1] = singleQuantileValues[1].first;
+                                        }
+                                    } else {
+                                        ++epochValues[1];
+                                    }
+                                }
+                                
+                                // Translate the result points to the 'original' domain
+                                for (auto& p : resultPoints) {
+                                    std::vector<ValueType> convertedP;
+                                    for (uint64_t i = 0; i < 2; ++i) {
+                                        if (p[i] == infinity) {
+                                            convertedP.push_back(storm::utility::infinity<ValueType>());
+                                        } else {
+                                            if (lowerCostBounds[i]) {
+                                                // Translate > to >=
+                                                ++p[i];
+                                            }
+                                            if (i == 1 && storm::solver::maximize(optimizationDirections[i])) {
+                                                // When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction
+                                                --p[i];
+                                            }
+                                            if (p[i] < 0) {
+                                                // Skip this point
+                                                convertedP.clear();
+                                                continue;
+                                            }
+                                            convertedP.push_back(storm::utility::convertNumber<ValueType>(p[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor);
+                                        }
+                                    }
+                                    if (!convertedP.empty()) {
+                                        result.push_back(convertedP);
+                                    }
+                                }
+                                filterDominatedPoints(result, optimizationDirections);
+                            } else if (maxmaxProbSatisfiesFormula) {
+                                // i.e., all bound values satisfy the formula
+                                if (storm::solver::minimize(optimizationDirections[0])) {
+                                    result = {{storm::utility::zero<ValueType>(), storm::utility::zero<ValueType>()}};
+                                } else {
+                                    result = {{ storm::utility::infinity<ValueType>(),  storm::utility::infinity<ValueType>()}};
+                                }
+                            } else {
+                                // i.e., no bound value satisfies the formula
+                                result = {{}};
+                            }
+                        } else {
+                            // TODO: this is an "unreasonable" case
+                        }
+                    } else {
+                        // TODO: find reasonable min/max cases
+                    }
+                    return result;
+                }
+                
                 template<typename ModelType>
                 std::pair<uint64_t, typename ModelType::ValueType> QuantileHelper<ModelType>::computeQuantileForDimension(Environment const& env, uint64_t dimension) const {
-                    // We assume that their is one bound value that violates the quantile and one bound value that satisfies it.
+                    // We assume that there is one bound value that violates the quantile and one bound value that satisfies it.
                     
-                    // 'Drop' all other open bounds
+                    // Let all other open bounds approach infinity
                     std::vector<BoundTransformation> bts(getDimension(), BoundTransformation::None);
+                    storm::storage::BitVector otherLowerBoundedDimensions(getDimension(), false);
                     for (auto const& d : getOpenDimensions()) {
                         if (d != dimension) {
-                            bts[d] = BoundTransformation::GreaterEqualZero;
+                            if (quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasLowerBound(d)) {
+                                bts[d] = BoundTransformation::GreaterZero;
+                                otherLowerBoundedDimensions.set(d, true);
+                            } else {
+                                bts[d] = BoundTransformation::GreaterEqualZero;
+                            }
                         }
                     }
+                    
+                    bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
+                    bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
+                    
+                    storm::storage::BitVector nonMecChoices;
+                    if (!otherLowerBoundedDimensions.empty()) {
+                        STORM_LOG_ASSERT(!upperCostBound, "It is not possible to let other open dimensions approach infinity if this is an upper cost bound and others are lower cost bounds.");
+                        // Get the choices that do not lie on a mec
+                        nonMecChoices.resize(model.getNumberOfChoices(), true);
+                        auto mecDecomposition = storm::storage::MaximalEndComponentDecomposition<ValueType>(model);
+                        for (auto const& mec : mecDecomposition) {
+                            for (auto const& stateChoicesPair : mec) {
+                                for (auto const& choice : stateChoicesPair.second) {
+                                    nonMecChoices.set(choice, false);
+                                }
+                            }
+                        }
+                    }
+                    
                     auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts);
-                    MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula);
+                    MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, [&](std::vector<EpochManager::Epoch>& epochSteps, EpochManager const& epochManager) {
+                        if (!otherLowerBoundedDimensions.empty()) {
+                            for (auto const& choice : nonMecChoices) {
+                                for (auto const& dim : otherLowerBoundedDimensions) {
+                                    epochManager.setDimensionOfEpoch(epochSteps[choice], dim, 0);
+                                }
+                            }
+                        }
+                    });
 
-                    bool upperCostBound = transformedFormula->getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
-                   // bool lowerProbabilityBound = storm::logic::isLowerBound(transformedFormula->getBound().comparisonType);
-                    bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
                     
                     // initialize data that will be needed for each epoch
                     auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
@@ -206,7 +409,7 @@ namespace storm {
                     std::vector<ValueType> x, b;
                     std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
                     
-                    int64_t epochValue = 0; // -1 is actally allowed since we express >=0 as >-1
+                    uint64_t epochValue = 0;
                     std::set<typename EpochManager::Epoch> exploredEpochs;
                     while (true) {
                         auto currEpoch = rewardUnfolding.getStartEpoch(true);
@@ -225,10 +428,20 @@ namespace storm {
                         // If the reward bound should be as small as possible, we should stop as soon as the property is satisfied.
                         // If the reward bound should be as large as possible, we should stop as soon as sthe property is violated and then go a step backwards
                         if (minimizingRewardBound && propertySatisfied) {
+                            // We found a satisfying value!
+                            if (!upperCostBound) {
+                                // The rewardunfolding assumes strict lower cost bounds while we assume non-strict ones. Hence, >B becomes >=(B+1).
+                                ++epochValue;
+                            }
                             break;
                         } else if (!minimizingRewardBound && !propertySatisfied) {
-                            STORM_LOG_ASSERT(epochValue > 0 || !upperCostBound, "The property does not seem to be satisfiable. This case should have been treated earlier.");
-                            --epochValue;
+                            // We found a non-satisfying value. Go one step back to get the largest satisfying value.
+                            // ... however, lower cost bounds need to be converted from strict  bounds >B to non strict bounds >=(B+1).
+                            // Hence, no operation is necessary in this case.
+                            if (upperCostBound) {
+                                STORM_LOG_ASSERT(epochValue > 0, "The property does not seem to be satisfiable. This case should have been excluded earlier.");
+                                --epochValue;
+                            }
                             break;
                         }
                         ++epochValue;
@@ -278,6 +491,7 @@ namespace storm {
                             }
                         }
                     }
+                    std::cout << "nonmec choices are " << nonMecChoices << std::endl;
 
                     MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, [&](std::vector<EpochManager::Epoch>& epochSteps, EpochManager const& epochManager) {
                         if (!minimizingLowerBoundedDimensions.empty()) {
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index 789e1d727..a0cd22fdb 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -22,6 +22,7 @@ namespace storm {
                     std::vector<std::vector<ValueType>> computeQuantile(Environment const& env);
 
                 private:
+                    std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment const& env);
                     
                     /*!
                      * Computes the infimum over bound values in minimizing dimensions / the supremum over bound values in maximizing dimensions