From 8c31019bd559f1625b9e24b745168eb1260dc140 Mon Sep 17 00:00:00 2001 From: TimQu Date: Sun, 13 Nov 2016 19:03:22 +0100 Subject: [PATCH] improved selection of digitization constant Former-commit-id: 72306b1879c92b5a5b5ca5b35e34161721a0175a --- .../pcaa/SparseMaPcaaWeightVectorChecker.cpp | 29 ++++----- .../pcaa/SparseMaPcaaWeightVectorChecker.h | 4 +- .../pcaa/SparsePcaaAchievabilityQuery.cpp | 47 ++++++++------ .../pcaa/SparsePcaaAchievabilityQuery.h | 5 ++ .../pcaa/SparsePcaaParetoQuery.cpp | 13 ++-- .../pcaa/SparsePcaaQuantitativeQuery.cpp | 63 +++++++++++++------ .../pcaa/SparsePcaaQuantitativeQuery.h | 6 ++ .../multiobjective/pcaa/SparsePcaaQuery.cpp | 27 ++------ .../pcaa/SparsePcaaWeightVectorChecker.cpp | 14 +++-- .../pcaa/SparsePcaaWeightVectorChecker.h | 41 +++++++----- src/storage/geometry/Halfspace.h | 28 +++++++-- src/storage/geometry/Polytope.cpp | 2 +- 12 files changed, 169 insertions(+), 110 deletions(-) diff --git a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp index 324cd1396..09806882c 100644 --- a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp +++ b/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 template ::SupportsExponential, int>::type> - VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant() const { + VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant(std::vector 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 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(); for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { auto const& obj = this->objectives[objIndex]; VT precisionOfObj = storm::utility::zero(); @@ -201,11 +205,9 @@ namespace storm { if(obj.upperTimeBound) { precisionOfObj += storm::utility::one() - (eToPowerOfMinusMaxRateTimesBound[objIndex].second * storm::utility::pow(storm::utility::one() + 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 template ::SupportsExponential, int>::type> - VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant() const { + VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant(std::vector 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() + 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>; - template double SparseMaPcaaWeightVectorChecker>::getDigitizationConstant() const; + template double SparseMaPcaaWeightVectorChecker>::getDigitizationConstant(std::vector const& direction) const; template void SparseMaPcaaWeightVectorChecker>::digitize(SparseMaPcaaWeightVectorChecker>::SubModel& subModel, double const& digitizationConstant) const; template void SparseMaPcaaWeightVectorChecker>::digitizeTimeBounds(SparseMaPcaaWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant); #ifdef STORM_HAVE_CARL // template class SparseMaPcaaWeightVectorChecker>; - // template storm::RationalNumber SparseMaPcaaWeightVectorChecker>::getDigitizationConstant() const; + // template storm::RationalNumber SparseMaPcaaWeightVectorChecker>::getDigitizationConstant(std::vector const& direction) const; // template void SparseMaPcaaWeightVectorChecker>::digitize(SparseMaPcaaWeightVectorChecker>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const; // template void SparseMaPcaaWeightVectorChecker>::digitizeTimeBounds(SparseMaPcaaWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant); #endif diff --git a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h index 24166e023..0540b39cc 100644 --- a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h +++ b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h @@ -92,9 +92,9 @@ namespace storm { * Retrieves the delta used for digitization */ template ::SupportsExponential, int>::type = 0> - VT getDigitizationConstant() const; + VT getDigitizationConstant(std::vector const& weightVector) const; template ::SupportsExponential, int>::type = 0> - VT getDigitizationConstant() const; + VT getDigitizationConstant(std::vector const& weightVector) const; /*! * Digitizes the given matrix and vectors w.r.t. the given digitization constant and the given rate vector. diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp index a330208db..fa0791afb 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp @@ -21,14 +21,8 @@ namespace storm { STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType::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(storm::settings::getModule().getPrecision()); - gap /= (storm::utility::one() + storm::utility::one()); - gap /= storm::utility::sqrt(static_cast(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(0.1)); } template @@ -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 + void SparsePcaaAchievabilityQuery::updateWeightedPrecision(WeightVector const& weights) { + // Our heuristic considers the distance between the under- and the over approximation w.r.t. the given direction + std::pair optimizationResOverApprox = this->overApproximation->optimize(weights); + if(optimizationResOverApprox.second) { + std::pair 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(), "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(distance)); + } + } + // do not update the precision if one of the approximations is unbounded in the provided direction + } template bool SparsePcaaAchievabilityQuery::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope) { std::vector> halfspaces = polytope->getHalfspaces(); for(auto const& h : halfspaces) { - GeometryValueType distance = h.distance(thresholds); - if(distance < storm::utility::zero()) { - return false; - } - if(distance == storm::utility::zero()) { - // 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()) { - 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()) { + return false; + } } } + } else { + return false; } } return true; diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h index 0b4ca9417..638e2174c 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h +++ b/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. diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.cpp index 14605ad28..e8746441f 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.cpp +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.cpp @@ -20,13 +20,10 @@ namespace storm { SparsePcaaParetoQuery::SparsePcaaParetoQuery(SparsePcaaPreprocessorReturnType& preprocessorResult) : SparsePcaaQuery(preprocessorResult) { STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType::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(storm::settings::getModule().getPrecision()); - gap /= (storm::utility::one() + storm::utility::one()); - gap /= storm::utility::sqrt(static_cast(this->objectives.size())); - this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap); + // Set the precision of the weight vector checker + typename SparseModelType::ValueType weightedPrecision = storm::utility::convertNumber(storm::settings::getModule().getPrecision()); + weightedPrecision /= typename SparseModelType::ValueType(2); + this->weightVectorChecker->setWeightedPrecision(weightedPrecision); } @@ -68,7 +65,7 @@ namespace storm { GeometryValueType farestDistance = storm::utility::zero(); 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; diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp index d16fd9eb7..05b83b443 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp +++ b/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(storm::settings::getModule().getPrecision()); - gap /= (storm::utility::one() + storm::utility::one()); - gap /= storm::utility::sqrt(static_cast(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(0.1)); } template @@ -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 + void SparsePcaaQuantitativeQuery::updateWeightedPrecisionInAchievabilityPhase(WeightVector const& weights) { + // Our heuristic considers the distance between the under- and the over approximation w.r.t. the given direction + std::pair optimizationResOverApprox = this->overApproximation->optimize(weights); + if(optimizationResOverApprox.second) { + std::pair 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(), "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(distance)); + } + } + // do not update the precision if one of the approximations is unbounded in the provided direction + } template GeometryValueType SparsePcaaQuantitativeQuery::improveSolution() { @@ -130,28 +143,42 @@ namespace storm { thresholds[indexOfOptimizingObjective] = result + storm::utility::one(); } 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 + void SparsePcaaQuantitativeQuery::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(storm::settings::getModule().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(weightedPrecision)); + } + + template bool SparsePcaaQuantitativeQuery::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope) { std::vector> halfspaces = polytope->getHalfspaces(); for(auto const& h : halfspaces) { - GeometryValueType distance = h.distance(thresholds); - if(distance < storm::utility::zero()) { - return false; - } - if(distance == storm::utility::zero()) { - // 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()) { - 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()) { + return false; + } } } + } else { + return false; } } return true; diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h index b2313e529..b25cf5c18 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h +++ b/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 */ diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp index 93b17c126..f38c01669 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp +++ b/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> halfspaces = underApproximation->getHalfspaces(); uint_fast64_t farestHalfspaceIndex = halfspaces.size(); GeometryValueType farestDistance = -storm::utility::one(); 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 distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); + if(!storm::utility::isZero(distance)) { storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero(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::performRefinementStep(WeightVector&& direction) { // Normalize the direction vector so that the entries sum up to one storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero())); - // 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 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 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(distance) * storm::utility::convertNumber(0.5))); - break; - } - } weightVectorChecker->check(storm::utility::vector::convertNumericVector(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(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(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(h.invert().euclideanDistance(refinementSteps.back().upperBoundPoint)) << "."); } overApproximation = overApproximation->intersection(h); STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp index 2ee474dd9..9da82e9fa 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp +++ b/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(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(), 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 - void SparsePcaaWeightVectorChecker::setMaximumLowerUpperBoundGap(ValueType const& value) { - this->maximumLowerUpperBoundGap = value; + void SparsePcaaWeightVectorChecker::setWeightedPrecision(ValueType const& weightedPrecision) { + this->weightedPrecision = weightedPrecision; } template - typename SparsePcaaWeightVectorChecker::ValueType const& SparsePcaaWeightVectorChecker::getMaximumLowerUpperBoundGap() const { - return this->maximumLowerUpperBoundGap; + typename SparsePcaaWeightVectorChecker::ValueType const& SparsePcaaWeightVectorChecker::getWeightedPrecision() const { + return this->weightedPrecision; } template diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h index df69ec059..c13854818 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h @@ -47,16 +47,6 @@ namespace storm { */ void check(std::vector 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 getLowerBoundsOfInitialStateResults() const; std::vector 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> 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 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> 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 offsetsToLowerBound; std::vector offsetsToUpperBound; // The scheduler that maximizes the weighted rewards diff --git a/src/storage/geometry/Halfspace.h b/src/storage/geometry/Halfspace.h index 64c0cd8a0..2b5705049 100644 --- a/src/storage/geometry/Halfspace.h +++ b/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 const& point) const { - return offset() - storm::utility::vector::dotProduct(point, normalVector()); + return std::max(storm::utility::zero(), 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 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 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 invert() const { + std::vector resNormalVector = normalVector(); + storm::utility::vector::scaleVectorInPlace(resNormalVector, -storm::utility::one()); + return Halfspace(std::move(resNormalVector), -offset()); + } /* * Returns a string representation of this Halfspace. diff --git a/src/storage/geometry/Polytope.cpp b/src/storage/geometry/Polytope.cpp index 150535dd0..2228e00e8 100644 --- a/src/storage/geometry/Polytope.cpp +++ b/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