improved selection of digitization constant

Former-commit-id: 72306b1879

src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp

@@ -45,7 +45,7 @@ namespace storm {
                 SubModel PS = createSubModel(false, weightedRewardVector);
 
                 // Apply digitization to Markovian transitions
-                ValueType digitizationConstant = getDigitizationConstant();
+                ValueType digitizationConstant = getDigitizationConstant(weightVector);
                 digitize(MS, digitizationConstant);
                 
                 // Get for each occurring (digitized) timeBound the indices of the objectives with that bound.
@@ -145,13 +145,15 @@ namespace storm {
             
             template <class SparseMaModelType>
             template <typename VT, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type>
-            VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::getDigitizationConstant() const {
+            VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::getDigitizationConstant(std::vector<ValueType> const& weightVector) const {
                 STORM_LOG_DEBUG("Retrieving digitization constant");
                 // We need to find a delta such that for each objective it holds that lowerbound/delta , upperbound/delta are natural numbers and
-                // If there is a lower and an upper bound:
-                //     1 - e^(-maxRate lowerbound) * (1 + maxRate delta) ^ (lowerbound / delta) + 1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta) + (1-e^(-maxRate delta) <= maximumLowerUpperBoundGap
-                // If there is only an upper bound:
-                //     1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta) <= maximumLowerUpperBoundGap
+                // sum_{obj_i} (
+                //   If obj_i has a lower and an upper bound:
+                //     weightVector_i * (1 - e^(-maxRate lowerbound) * (1 + maxRate delta) ^ (lowerbound / delta) + 1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta) + (1-e^(-maxRate delta)))
+                //   If there is only an upper bound:
+                //     weightVector_i * ( 1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta))
+                // ) <= this->maximumLowerUpperDistance
                 
                 // Initialize some data for fast and easy access
                 VT const maxRate = this->model.getMaximalExitRate();
@@ -175,6 +177,7 @@ namespace storm {
                     // There are no time bounds. In this case, one is a valid digitization constant.
                     return storm::utility::one<VT>();
                 }
+                VT goalPrecisionTimesNorm = this->weightedPrecision * storm::utility::sqrt(storm::utility::vector::dotProduct(weightVector, weightVector));
                 
                 // We brute-force a delta, since a direct computation is apparently not easy.
                 // Also note that the number of times this loop runs is a lower bound for the number of minMaxSolver invocations.
@@ -191,6 +194,7 @@ namespace storm {
                         }
                     }
                     if(deltaValid) {
+                        VT weightedPrecisionForCurrentDelta = storm::utility::zero<VT>();
                         for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                             auto const& obj = this->objectives[objIndex];
                             VT precisionOfObj = storm::utility::zero<VT>();
@@ -201,11 +205,9 @@ namespace storm {
                             if(obj.upperTimeBound) {
                                 precisionOfObj += storm::utility::one<VT>() - (eToPowerOfMinusMaxRateTimesBound[objIndex].second * storm::utility::pow(storm::utility::one<VT>() + maxRate * delta, *obj.upperTimeBound / delta) );
                             }
-                            if(precisionOfObj > this->maximumLowerUpperBoundGap) {
-                                deltaValid = false;
-                                break;
-                            }
+                            weightedPrecisionForCurrentDelta += weightVector[objIndex] * precisionOfObj;
                         }
+                        deltaValid &= weightedPrecisionForCurrentDelta <= goalPrecisionTimesNorm;
                     }
                     if(deltaValid) {
                         break;
@@ -220,7 +222,7 @@ namespace storm {
             
             template <class SparseMaModelType>
             template <typename VT, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type>
-            VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::getDigitizationConstant() const {
+            VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::getDigitizationConstant(std::vector<ValueType> const& weightVector) const {
                   STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type.");
             }
             
@@ -282,7 +284,6 @@ namespace storm {
                         digitizationError -= std::exp(-maxRate * (*obj.upperTimeBound)) * storm::utility::pow(storm::utility::one<VT>() + maxRate * digitizationConstant, digitizedBound);
                         errorAwayFromZero += digitizationError;
                     }
-                    STORM_LOG_ASSERT(errorTowardsZero + errorAwayFromZero <= this->maximumLowerUpperBoundGap, "Precision not sufficient.");
                     if (obj.rewardsArePositive) {
                         this->offsetsToLowerBound[objIndex] = -errorTowardsZero;
                         this->offsetsToUpperBound[objIndex] = errorAwayFromZero;
@@ -408,12 +409,12 @@ namespace storm {
             
             
             template class SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>;
-            template double SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::getDigitizationConstant<double>() const;
+            template double SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::getDigitizationConstant<double>(std::vector<double> const& direction) const;
             template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitize<double>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::SubModel& subModel, double const& digitizationConstant) const;
             template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitizeTimeBounds<double>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant);
 #ifdef STORM_HAVE_CARL
 //            template class SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
-   //         template storm::RationalNumber SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::getDigitizationConstant<storm::RationalNumber>() const;
+   //         template storm::RationalNumber SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::getDigitizationConstant<storm::RationalNumber>(std::vector<double> const& direction) const;
  //           template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitize<storm::RationalNumber>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const;
 //            template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitizeTimeBounds<storm::RationalNumber>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant);
 #endif

src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h

@@ -92,9 +92,9 @@ namespace storm {
                  * Retrieves the delta used for digitization
                  */
                 template <typename VT = ValueType, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type = 0>
-                VT getDigitizationConstant() const;
+                VT getDigitizationConstant(std::vector<ValueType> const& weightVector) const;
                 template <typename VT = ValueType, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type = 0>
-                VT getDigitizationConstant() const;
+                VT getDigitizationConstant(std::vector<ValueType> const& weightVector) const;
                 
                 /*!
                  * Digitizes the given matrix and vectors w.r.t. the given digitization constant and the given rate vector.

src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp

@@ -21,14 +21,8 @@ namespace storm {
                 STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType<SparseModelType>::QueryType::Achievability, "Invalid query Type");
                 initializeThresholdData();
                 
-                // Set the maximum gap between lower and upper bound of the weightVectorChecker result.
-                // This is the maximal edge length of the box we have to consider around each computed point
-                // We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
-                typename SparseModelType::ValueType gap = storm::utility::convertNumber<typename SparseModelType::ValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
-                gap /= (storm::utility::one<typename SparseModelType::ValueType>() + storm::utility::one<typename SparseModelType::ValueType>());
-                gap /= storm::utility::sqrt(static_cast<typename SparseModelType::ValueType>(this->objectives.size()));
-                this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
-                
+                // Set the precision of the weight vector checker. Will be refined during the computation
+                this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(0.1));
             }
             
             template <class SparseModelType, typename GeometryValueType>
@@ -55,6 +49,7 @@ namespace storm {
                 // repeatedly refine the over/ under approximation until the threshold point is either in the under approx. or not in the over approx.
                 while(!this->maxStepsPerformed()){
                     WeightVector separatingVector = this->findSeparatingVector(thresholds);
+                    this->updateWeightedPrecision(separatingVector);
                     this->performRefinementStep(std::move(separatingVector));
                     if(!checkIfThresholdsAreSatisfied(this->overApproximation)){
                         return false;
@@ -67,23 +62,39 @@ namespace storm {
                 return false;
             }
 
+            template <class SparseModelType, typename GeometryValueType>
+            void SparsePcaaAchievabilityQuery<SparseModelType, GeometryValueType>::updateWeightedPrecision(WeightVector const& weights) {
+                // Our heuristic considers the distance between the under- and the over approximation w.r.t. the given direction
+                std::pair<Point, bool> optimizationResOverApprox = this->overApproximation->optimize(weights);
+                if(optimizationResOverApprox.second) {
+                    std::pair<Point, bool> optimizationResUnderApprox = this->underApproximation->optimize(weights);
+                    if(optimizationResUnderApprox.second) {
+                        GeometryValueType distance = storm::utility::vector::dotProduct(optimizationResOverApprox.first, weights) - storm::utility::vector::dotProduct(optimizationResUnderApprox.first, weights);
+                        STORM_LOG_ASSERT(distance >= storm::utility::zero<GeometryValueType>(), "Negative distance between under- and over approximation was not expected");
+                        // Normalize the distance by dividing it with the Euclidean Norm of the weight-vector
+                        distance /= storm::utility::sqrt(storm::utility::vector::dotProduct(weights, weights));
+                        distance /= GeometryValueType(2);
+                        this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(distance));
+                    }
+                }
+                // do not update the precision if one of the approximations is unbounded in the provided direction
+            }
             
             template <class SparseModelType, typename GeometryValueType>
             bool SparsePcaaAchievabilityQuery<SparseModelType, GeometryValueType>::checkIfThresholdsAreSatisfied(std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) {
                 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = polytope->getHalfspaces();
                 for(auto const& h : halfspaces) {
-                    GeometryValueType distance = h.distance(thresholds);
-                    if(distance < storm::utility::zero<GeometryValueType>()) {
-                        return false;
-                    }
-                    if(distance == storm::utility::zero<GeometryValueType>()) {
-                        // In this case, the thresholds point is on the boundary of the polytope.
-                        // Check if this is problematic for the strict thresholds
-                        for(auto strictThreshold : strictThresholds) {
-                            if(h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
-                                return false;
+                    if(storm::utility::isZero(h.distance(thresholds))) {
+                        // Check if the threshold point is on the boundary of the halfspace and whether this is violates strict thresholds
+                        if(h.isPointOnBoundary(thresholds)) {
+                            for(auto strictThreshold : strictThresholds) {
+                                if(h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
+                                    return false;
+                                }
                             }
                         }
+                    } else {
+                        return false;
                     }
                 }
                 return true;

src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h

@@ -40,6 +40,11 @@ namespace storm {
                  */
                 bool checkAchievability();
                 
+                /*
+                 * Updates the precision of the weightVectorChecker w.r.t. the provided weights
+                 */
+                void updateWeightedPrecision(WeightVector const& weights);
+                
                 /*
                  * Returns true iff there is one point in the given polytope that satisfies the given thresholds.
                  * It is assumed that the given polytope contains the downward closure of its vertices.

src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.cpp

@@ -20,13 +20,10 @@ namespace storm {
             SparsePcaaParetoQuery<SparseModelType, GeometryValueType>::SparsePcaaParetoQuery(SparsePcaaPreprocessorReturnType<SparseModelType>& preprocessorResult) : SparsePcaaQuery<SparseModelType, GeometryValueType>(preprocessorResult) {
                 STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType<SparseModelType>::QueryType::Pareto, "Invalid query Type");
                 
-                // Set the maximum gap between lower and upper bound of the weightVectorChecker result.
-                // This is the maximal edge length of the box we have to consider around each computed point
-                // We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
-                typename SparseModelType::ValueType gap = storm::utility::convertNumber<typename SparseModelType::ValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
-                gap /= (storm::utility::one<typename SparseModelType::ValueType>() + storm::utility::one<typename SparseModelType::ValueType>());
-                gap /= storm::utility::sqrt(static_cast<typename SparseModelType::ValueType>(this->objectives.size()));
-                this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
+                // Set the precision of the weight vector checker
+                typename SparseModelType::ValueType weightedPrecision = storm::utility::convertNumber<typename SparseModelType::ValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
+                weightedPrecision /= typename SparseModelType::ValueType(2);
+                this->weightVectorChecker->setWeightedPrecision(weightedPrecision);
                 
             }
             
@@ -68,7 +65,7 @@ namespace storm {
                     GeometryValueType farestDistance = storm::utility::zero<GeometryValueType>();
                     for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < underApproxHalfspaces.size(); ++halfspaceIndex) {
                         for(auto const& vertex : overApproxVertices) {
-                            GeometryValueType distance = -underApproxHalfspaces[halfspaceIndex].euclideanDistance(vertex);
+                            GeometryValueType distance = underApproxHalfspaces[halfspaceIndex].euclideanDistance(vertex);
                             if(distance > farestDistance) {
                                 farestHalfspaceIndex = halfspaceIndex;
                                 farestDistance = distance;

src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp

@@ -24,14 +24,8 @@ namespace storm {
                 indexOfOptimizingObjective = *preprocessorResult.indexOfOptimizingObjective;
                 initializeThresholdData();
                 
-                // Set the maximum gap between lower and upper bound of the weightVectorChecker result.
-                // This is the maximal edge length of the box we have to consider around each computed point
-                // We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
-                typename SparseModelType::ValueType gap = storm::utility::convertNumber<typename SparseModelType::ValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
-                gap /= (storm::utility::one<typename SparseModelType::ValueType>() + storm::utility::one<typename SparseModelType::ValueType>());
-                gap /= storm::utility::sqrt(static_cast<typename SparseModelType::ValueType>(this->objectives.size()));
-                this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
-                
+                // Set the maximum distance between lower and upper bound of the weightVectorChecker result.
+                this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(0.1));
             }
             
             template <class SparseModelType, typename GeometryValueType>
@@ -84,6 +78,7 @@ namespace storm {
                 
                 while(!this->maxStepsPerformed()){
                     WeightVector separatingVector = this->findSeparatingVector(thresholds);
+                    this->updateWeightedPrecisionInAchievabilityPhase(separatingVector);
                     this->performRefinementStep(std::move(separatingVector));
                     //Pick the threshold for the optimizing objective low enough so valid solutions are not excluded
                     thresholds[indexOfOptimizingObjective] = std::min(thresholds[indexOfOptimizingObjective], this->refinementSteps.back().lowerBoundPoint[indexOfOptimizingObjective]);
@@ -97,6 +92,24 @@ namespace storm {
                 STORM_LOG_ERROR("Could not check whether thresholds are achievable: Exceeded maximum number of refinement steps");
                 return false;
             }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInAchievabilityPhase(WeightVector const& weights) {
+                // Our heuristic considers the distance between the under- and the over approximation w.r.t. the given direction
+                std::pair<Point, bool> optimizationResOverApprox = this->overApproximation->optimize(weights);
+                if(optimizationResOverApprox.second) {
+                    std::pair<Point, bool> optimizationResUnderApprox = this->underApproximation->optimize(weights);
+                    if(optimizationResUnderApprox.second) {
+                        GeometryValueType distance = storm::utility::vector::dotProduct(optimizationResOverApprox.first, weights) - storm::utility::vector::dotProduct(optimizationResUnderApprox.first, weights);
+                        STORM_LOG_ASSERT(distance >= storm::utility::zero<GeometryValueType>(), "Negative distance between under- and over approximation was not expected");
+                        // Normalize the distance by dividing it with the Euclidean Norm of the weight-vector
+                        distance /= storm::utility::sqrt(storm::utility::vector::dotProduct(weights, weights));
+                        distance /= GeometryValueType(2);
+                        this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(distance));
+                    }
+                }
+                // do not update the precision if one of the approximations is unbounded in the provided direction
+            }
 
             template <class SparseModelType, typename GeometryValueType>
             GeometryValueType SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::improveSolution() {
@@ -130,28 +143,42 @@ namespace storm {
                         thresholds[indexOfOptimizingObjective] = result + storm::utility::one<GeometryValueType>();
                     }
                     WeightVector separatingVector = this->findSeparatingVector(thresholds);
+                    this->updateWeightedPrecisionInImprovingPhase(separatingVector);
                     this->performRefinementStep(std::move(separatingVector));
                 }
                STORM_LOG_ERROR("Could not reach the desired precision: Exceeded maximum number of refinement steps");
                return result;
             }
             
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInImprovingPhase(WeightVector const& weights) {
+                STORM_LOG_THROW(!storm::utility::isZero(weights[this->indexOfOptimizingObjective]), exceptions::UnexpectedException, "The chosen weight-vector gives zero weight for the objective that is to be optimized.");
+                // If weighs[indexOfOptimizingObjective] is low, the computation of the weightVectorChecker needs to be more precise.
+                // Our heuristic ensures that if p is the new vertex of the under-approximation, then max{ eps | p' = p + (0..0 eps 0..0) is in the over-approximation } <= multiobjective_precision/2
+                    GeometryValueType weightedPrecision = weights[this->indexOfOptimizingObjective] * storm::utility::convertNumber<GeometryValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
+                    // Normalize by division with the Euclidean Norm of the weight-vector
+                    weightedPrecision /= storm::utility::sqrt(storm::utility::vector::dotProduct(weights, weights));
+                    weightedPrecision /= GeometryValueType(2);
+                    this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(weightedPrecision));
+            }
+            
+            
             template <class SparseModelType, typename GeometryValueType>
             bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::checkIfThresholdsAreSatisfied(std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) {
                 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = polytope->getHalfspaces();
                 for(auto const& h : halfspaces) {
-                    GeometryValueType distance = h.distance(thresholds);
-                    if(distance < storm::utility::zero<GeometryValueType>()) {
-                        return false;
-                    }
-                    if(distance == storm::utility::zero<GeometryValueType>()) {
-                        // In this case, the thresholds point is on the boundary of the polytope.
-                        // Check if this is problematic for the strict thresholds
-                        for(auto strictThreshold : strictThresholds) {
-                            if(h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
-                                return false;
+                    if(storm::utility::isZero(h.distance(thresholds))) {
+                        // Check if the threshold point is on the boundary of the halfspace and whether this is violates strict thresholds
+                        if(h.isPointOnBoundary(thresholds)) {
+                            for(auto strictThreshold : strictThresholds) {
+                                if(h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
+                                    return false;
+                                }
                             }
                         }
+                    } else {
+                        return false;
                     }
                 }
                 return true;

src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h

@@ -40,6 +40,12 @@ namespace storm {
                  */
                 bool checkAchievability();
                 
+                /*
+                 * Updates the precision of the weightVectorChecker w.r.t. the provided weights
+                 */
+                void updateWeightedPrecisionInAchievabilityPhase(WeightVector const& weights);
+                void updateWeightedPrecisionInImprovingPhase(WeightVector const& weights);
+                
                 /*
                  * Given that the thresholds are achievable, this function further refines the approximations and returns the optimized value
                  */

src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp

@@ -61,13 +61,14 @@ namespace storm {
                 
                 // Reaching this point means that the underApproximation contains halfspaces. The seperating vector has to be the normal vector of one of these halfspaces.
                 // We pick one with maximal distance to the given point. However, Dirac weight vectors that only assign a non-zero weight to a single objective take precedence.
+                STORM_LOG_ASSERT(!underApproximation->contains(pointToBeSeparated), "Tried to find a separating point but the point is already contained in the underApproximation");
                 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = underApproximation->getHalfspaces();
                 uint_fast64_t farestHalfspaceIndex = halfspaces.size();
                 GeometryValueType farestDistance = -storm::utility::one<GeometryValueType>();
                 bool foundSeparatingDiracVector = false;
                 for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) {
-                    GeometryValueType distance = -halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated);
-                    if(distance >= storm::utility::zero<GeometryValueType>()) {
+                    GeometryValueType distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated);
+                    if(!storm::utility::isZero(distance)) {
                         storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero<GeometryValueType>(halfspaces[halfspaceIndex].normalVector());
                         bool isSingleObjectiveVector = nonZeroVectorEntries.getNumberOfSetBits() == 1 && diracWeightVectorsToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0));
                         // Check if this halfspace is a better candidate than the current one
@@ -91,27 +92,7 @@ namespace storm {
             void SparsePcaaQuery<SparseModelType, GeometryValueType>::performRefinementStep(WeightVector&& direction) {
                 // Normalize the direction vector so that the entries sum up to one
                 storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one<GeometryValueType>() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero<GeometryValueType>()));
-                // Check if we already did a refinement step with that direction vector. If this is the case, we increase the precision.
-                // We start with the most recent steps to consider the most recent result for this direction vector
-                boost::optional<typename SparseModelType::ValueType> oldMaximumLowerUpperBoundGap;
-                for(auto stepIt = refinementSteps.rbegin(); stepIt != refinementSteps.rend(); ++stepIt) {
-                    if(stepIt->weightVector == direction) {
-                        STORM_LOG_WARN("Performing multiple refinement steps with the same direction vector.");
-                        oldMaximumLowerUpperBoundGap = weightVectorChecker->getMaximumLowerUpperBoundGap();
-                        std::vector<GeometryValueType> lowerUpperDistances = stepIt->upperBoundPoint;
-                        storm::utility::vector::subtractVectors(lowerUpperDistances, stepIt->lowerBoundPoint, lowerUpperDistances);
-                        // shorten the distance between lower and upper bound for the new result by multiplying the current distance with 0.5
-                        // TODO: try other values/strategies?
-                        GeometryValueType distance = storm::utility::sqrt(storm::utility::vector::dotProduct(lowerUpperDistances, lowerUpperDistances));
-                        weightVectorChecker->setMaximumLowerUpperBoundGap(std::min(*oldMaximumLowerUpperBoundGap, storm::utility::convertNumber<typename SparseModelType::ValueType>(distance) * storm::utility::convertNumber<typename SparseModelType::ValueType>(0.5)));
-                        break;
-                    }
-                }
                 weightVectorChecker->check(storm::utility::vector::convertNumericVector<typename SparseModelType::ValueType>(direction));
-                if(oldMaximumLowerUpperBoundGap) {
-                    // Reset the precision back to the previous values
-                    weightVectorChecker->setMaximumLowerUpperBoundGap(*oldMaximumLowerUpperBoundGap);
-                }
                 STORM_LOG_DEBUG("weighted objectives checker result (lower bounds) is " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(weightVectorChecker->getLowerBoundsOfInitialStateResults())));
                 RefinementStep step;
                 step.weightVector = direction;
@@ -136,7 +117,7 @@ namespace storm {
                 if(maximumOffset > h.offset()){
                     // We correct the issue by shifting the halfspace such that it contains the underapproximation
                     h.offset() = maximumOffset;
-                    STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by " << storm::utility::convertNumber<double>(h.euclideanDistance(refinementSteps.back().upperBoundPoint)) << ".");
+                    STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by " << storm::utility::convertNumber<double>(h.invert().euclideanDistance(refinementSteps.back().upperBoundPoint)) << ".");
                 }
                 overApproximation = overApproximation->intersection(h);
                 STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true));

src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp

@@ -14,6 +14,7 @@
 #include "src/utility/vector.h"
 
 #include "src/exceptions/IllegalFunctionCallException.h"
+#include "src/exceptions/UnexpectedException.h"
 #include "src/exceptions/NotImplementedException.h"
 
 namespace storm {
@@ -64,16 +65,21 @@ namespace storm {
                     }
                 }
                 STORM_LOG_DEBUG("Weight vector check done. Lower bounds for results in initial state: " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(getLowerBoundsOfInitialStateResults())));
+                // Validate that the results are sufficiently precise
+                ValueType resultingWeightedPrecision = storm::utility::vector::dotProduct(getUpperBoundsOfInitialStateResults(), weightVector) - storm::utility::vector::dotProduct(getLowerBoundsOfInitialStateResults(), weightVector);
+                STORM_LOG_THROW(resultingWeightedPrecision >= storm::utility::zero<ValueType>(), storm::exceptions::UnexpectedException, "The distance between the lower and the upper result is negative.");
+                resultingWeightedPrecision /= storm::utility::sqrt(storm::utility::vector::dotProduct(weightVector, weightVector));
+                STORM_LOG_THROW(resultingWeightedPrecision <= weightedPrecision, storm::exceptions::UnexpectedException, "The desired precision was not reached");
             }
             
             template <class SparseModelType>
-            void SparsePcaaWeightVectorChecker<SparseModelType>::setMaximumLowerUpperBoundGap(ValueType const& value) {
-                this->maximumLowerUpperBoundGap = value;
+            void SparsePcaaWeightVectorChecker<SparseModelType>::setWeightedPrecision(ValueType const& weightedPrecision) {
+                this->weightedPrecision = weightedPrecision;
             }
             
             template <class SparseModelType>
-            typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType const& SparsePcaaWeightVectorChecker<SparseModelType>::getMaximumLowerUpperBoundGap() const {
-                return this->maximumLowerUpperBoundGap;
+            typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType const& SparsePcaaWeightVectorChecker<SparseModelType>::getWeightedPrecision() const {
+                return this->weightedPrecision;
             }
             
             template <class SparseModelType>

src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h

@@ -47,16 +47,6 @@ namespace storm {
                  */
                 void check(std::vector<ValueType> const& weightVector);
                 
-                /*!
-                 * Sets the maximum gap that is allowed between the lower and upper bound of the result of some objective.
-                 */
-                void setMaximumLowerUpperBoundGap(ValueType const& value);
-                
-                /*!
-                 * Retrieves the maximum gap that is allowed between the lower and upper bound of the result of some objective.
-                 */
-                ValueType const& getMaximumLowerUpperBoundGap() const;
-                
                 /*!
                  * 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.
@@ -66,6 +56,22 @@ namespace storm {
                 std::vector<ValueType> getLowerBoundsOfInitialStateResults() const;
                 std::vector<ValueType> getUpperBoundsOfInitialStateResults() const;
                 
+                
+                /*!
+                 * Sets the precision of this weight vector checker. After calling check() the following will hold:
+                 * Let h_lower and h_upper be two hyperplanes such that
+                 * * the normal vector is the provided weight-vector
+                 * * getLowerBoundsOfInitialStateResults() lies on h_lower and
+                 * * getUpperBoundsOfInitialStateResults() lies on h_upper.
+                 * Then the distance between the two hyperplanes is at most weightedPrecision
+                 */
+                void setWeightedPrecision(ValueType const& weightedPrecision);
+                
+                /*!
+                 * Returns the precision of this weight vector checker.
+                 */
+                ValueType const& getWeightedPrecision() const;
+                
                 /*!
                  * Retrieves a scheduler that induces the current values
                  * Note that check(..) has to be called before retrieving the scheduler. Otherwise, an exception is thrown.
@@ -128,19 +134,22 @@ namespace storm {
                 storm::storage::BitVector objectivesWithNoUpperTimeBound;
                 // stores the (discretized) state action rewards for each objective.
                 std::vector<std::vector<ValueType>> discreteActionRewards;
-                // stores the maximum gap that is allowed between the lower and upper bound of the result of some objective.
-                ValueType maximumLowerUpperBoundGap;
-
+                /* stores the precision of this weight vector checker. After calling check() the following will hold:
+                * Let h_lower and h_upper be two hyperplanes such that
+                * * the normal vector is the provided weight-vector
+                * * getLowerBoundsOfInitialStateResults() lies on h_lower and
+                * * getUpperBoundsOfInitialStateResults() lies on h_upper.
+                * Then the distance between the two hyperplanes is at most weightedPrecision */
+                ValueType weightedPrecision;
                 // Memory for the solution of the most recent call of check(..)
                 // becomes true after the first call of check(..)
                 bool checkHasBeenCalled;
                 // The result for the weighted reward vector (for all states of the model)
                 std::vector<ValueType> weightedResult;
-                // The lower bounds of the results for the individual objectives (w.r.t. all states of the model)
+                // The results for the individual objectives (w.r.t. all states of the model)
                 std::vector<std::vector<ValueType>> objectiveResults;
                 // Stores for each objective the distance between the computed result (w.r.t. the initial state) and a lower/upper bound for the actual result.
-                // The distances are stored as a (possibly negative) offset that has to be added to to the objectiveResults.
-                // Note that there is no guarantee that the lower/upper bounds are sound as long as the underlying solution method is not sound (e.g. standard value iteration).
+                // The distances are stored as a (possibly negative) offset that has to be added (+) to to the objectiveResults.
                 std::vector<ValueType> offsetsToLowerBound;
                 std::vector<ValueType> offsetsToUpperBound;
                 // The scheduler that maximizes the weighted rewards

src/storage/geometry/Halfspace.h

@@ -35,24 +35,40 @@ namespace storm {
                 }
                 
                 /*
-                 * Returns the (scaled) distance of the given point from the hyperplane defined by normalVector*x = offset.
-                 * The point is inside this halfspace iff the returned value is >=0
+                 * Returns the (scaled) distance of the given point from this halfspace.
+                 * If the point is inside this halfspace, the distance is 0.
                  * The returned value is the euclidean distance times the 2-norm of the normalVector.
                  * In contrast to the euclideanDistance method, there are no inaccuracies introduced (providing ValueType is exact for +, -, and *)
                  */
                 ValueType distance(std::vector<ValueType> const& point) const {
-                    return offset() - storm::utility::vector::dotProduct(point, normalVector());
+                    return std::max(storm::utility::zero<ValueType>(), storm::utility::vector::dotProduct(point, normalVector()) - offset());
                 }
                 
                 /*
-                 * Returns the euclidean distance of the point from the hyperplane defined by this.
-                 * The point is inside this halfspace iff the distance is >=0
+                 * Returns the euclidean distance of the point from this halfspace.
+                 * If the point is inside this halfspace, the distance is 0.
                  * Note that the euclidean distance is in general not a rational number (which can introduce inaccuracies).
                  */
                 ValueType euclideanDistance(std::vector<ValueType> const& point) const {
-                    // divide with the 2-norm of the normal vector
+                    // divide the distance with the 2-norm of the normal vector
                     return distance(point) / storm::utility::sqrt(storm::utility::vector::dotProduct(normalVector(), normalVector()));
                 }
+                
+                /*
+                 * Returns true iff the given point lies on the boundary of this halfspace (i.e., on the hyperplane given by normalVector()*x =offset
+                 */
+                bool isPointOnBoundary(std::vector<ValueType> const& point) const {
+                    return storm::utility::vector::dotProduct(point, normalVector()) == offset();
+                }
+                
+                /*
+                 * Returns the inverted Halfspace of this which represents the set (R^n \ this) union { x | x is on the boundary of this}
+                 */
+                Halfspace<ValueType> invert() const {
+                    std::vector<ValueType> resNormalVector = normalVector();
+                    storm::utility::vector::scaleVectorInPlace(resNormalVector, -storm::utility::one<ValueType>());
+                    return Halfspace<ValueType>(std::move(resNormalVector), -offset());
+                }
 
                 /*
                  * Returns a string representation of this Halfspace.

src/storage/geometry/Polytope.cpp

@@ -112,7 +112,7 @@ namespace storm {
                 for(auto const& h : halfspaces) {
                     storm::storage::BitVector verticesOnHalfspace(vertices.size(), false);
                     for(uint_fast64_t v = 0; v<vertices.size(); ++v) {
-                        if(storm::utility::isZero(h.distance(vertices[v]))) {
+                        if(h.isPointOnBoundary(vertices[v])) {
                             verticesOnHalfspace.set(v);
                         }
                     }