From 2bab103a87d3321bb3703390303724735f58956c Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 8 Jun 2016 12:41:47 +0200 Subject: [PATCH] numerical and pareto queries Former-commit-id: 7846df4ebee82aa874de07fb68e0677f0aa4837c --- .../SparseMdpMultiObjectiveModelChecker.cpp | 3 +- ...seMdpMultiObjectivePreprocessingHelper.cpp | 2 + ...SparseMultiObjectiveModelCheckerHelper.cpp | 353 ++++++++++++++++-- .../SparseMultiObjectiveModelCheckerHelper.h | 18 +- ...rseMultiObjectiveModelCheckerInformation.h | 4 +- src/storage/geometry/Halfspace.h | 22 +- 6 files changed, 359 insertions(+), 43 deletions(-) diff --git a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp index 145b3b4d8..d4f7920a9 100644 --- a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp +++ b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp @@ -41,8 +41,7 @@ namespace storm { info.printInformationToStream(std::cout); #ifdef STORM_HAVE_CARL - helper::SparseMultiObjectiveModelCheckerHelper helper(info); - helper.achievabilityQuery(); + helper::SparseMultiObjectiveModelCheckerHelper::check(info); #else STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi objective model checking currently requires carl."); #endif diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp index c5b44cf2b..9190b3608 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp @@ -61,9 +61,11 @@ namespace storm { storm::logic::OperatorFormula const& opFormula = formula->asOperatorFormula(); if(opFormula.hasBound()){ objective.threshold = opFormula.getBound().threshold; + objective.thresholdIsStrict = storm::logic::isStrict(opFormula.getBound().comparisonType); // Note that we minimize if the comparison type is an upper bound since we are interested in the EXISTENCE of a scheduler. objective.originalFormulaMinimizes = !storm::logic::isLowerBound(opFormula.getBound().comparisonType); } else if (opFormula.hasOptimalityType()){ + objective.thresholdIsStrict = false; objective.originalFormulaMinimizes = storm::solver::minimize(opFormula.getOptimalityType()); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Current objective " << opFormula << " does not specify whether to minimize or maximize"); diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp index 3af2715e3..c57b41215 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp @@ -12,6 +12,27 @@ namespace storm { namespace modelchecker { namespace helper { + template + void SparseMultiObjectiveModelCheckerHelper::check(Information & info) { + SparseMultiObjectiveModelCheckerHelper helper(info); + uint_fast64_t numOfObjectivesWithoutThreshold = 0; + for(auto const& obj : info.objectives) { + if(!obj.threshold) { + ++numOfObjectivesWithoutThreshold; + } + } + + if(numOfObjectivesWithoutThreshold == 0) { + helper.achievabilityQuery(); + } else if (numOfObjectivesWithoutThreshold == 1) { + helper.numericalQuery(); + } else if (numOfObjectivesWithoutThreshold == info.objectives.size()) { + helper.paretoQuery(); + } else { + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The number of objecties without threshold is not valid. It should be either 0 (achievabilityQuery), 1 (numericalQuery), or " << info.objectives.size() << " (paretoQuery). Got " << numOfObjectivesWithoutThreshold << " instead."); + } + } + template SparseMultiObjectiveModelCheckerHelper::SparseMultiObjectiveModelCheckerHelper(Information& info) : info(info), weightedObjectivesChecker(info) { overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); @@ -25,27 +46,141 @@ namespace storm { template void SparseMultiObjectiveModelCheckerHelper::achievabilityQuery() { - Point queryPoint; - queryPoint.reserve(info.objectives.size()); - for(auto const& obj : info.objectives) { - STORM_LOG_THROW(obj.threshold, storm::exceptions::UnexpectedException, "Can not perform achievabilityQuery: No threshold given for at least one objective."); - queryPoint.push_back(storm::utility::convertNumber(*obj.threshold)); + Point thresholds; + thresholds.reserve(info.objectives.size()); + storm::storage::BitVector strictThresholds(info.objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + STORM_LOG_THROW(info.objectives[objIndex].threshold, storm::exceptions::UnexpectedException, "Can not perform achievabilityQuery: No threshold given for at least one objective."); + thresholds.push_back(storm::utility::convertNumber(*info.objectives[objIndex].threshold)); + strictThresholds.set(objIndex, info.objectives[objIndex].thresholdIsStrict); } + storm::storage::BitVector individualObjectivesToBeChecked(info.objectives.size(), true); while(true) { //TODO introduce convergence criterion? (like max num of iterations) - storm::storage::geometry::Halfspace separatingHalfspace = findSeparatingHalfspace(queryPoint, individualObjectivesToBeChecked); - refineSolution(separatingHalfspace.normalVector()); - if(!overApproximation->contains(queryPoint)){ + WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); + refineSolution(separatingVector); + if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ std::cout << "PROPERTY VIOLATED" << std::endl; return; } - if(underApproximation->contains(queryPoint)){ + if(checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)){ std::cout << "PROPERTY SATISFIED" << std::endl; return; } } } + template + void SparseMultiObjectiveModelCheckerHelper::numericalQuery() { + Point thresholds; + thresholds.reserve(info.objectives.size()); + storm::storage::BitVector strictThresholds(info.objectives.size()); + std::vector> thresholdConstraints; + thresholdConstraints.reserve(info.objectives.size()-1); + uint_fast64_t optimizingObjIndex = info.objectives.size(); + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + if(info.objectives[objIndex].threshold) { + thresholds.push_back(storm::utility::convertNumber(*info.objectives[objIndex].threshold)); + WeightVector normalVector(info.objectives.size(), storm::utility::zero()); + normalVector[objIndex] = -storm::utility::one(); + thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back()); + std::cout << "adding threshold constraint" << thresholdConstraints.back().toString(true) << std::endl; + strictThresholds.set(objIndex, info.objectives[objIndex].thresholdIsStrict); + } else { + STORM_LOG_ASSERT(optimizingObjIndex == info.objectives.size(), "Numerical Query invoked but there are multiple objectives without threshold"); + optimizingObjIndex = objIndex; + thresholds.push_back(storm::utility::zero()); + } + } + STORM_LOG_ASSERT(optimizingObjIndex < info.objectives.size(), "Numerical Query invoked but there are no objectives without threshold"); + auto thresholdsAsPolytope = storm::storage::geometry::Polytope::create(thresholdConstraints); + WeightVector directionOfOptimizingObjective(info.objectives.size()); + directionOfOptimizingObjective[optimizingObjIndex] = storm::utility::one(); + + // Try to find one valid solution + storm::storage::BitVector individualObjectivesToBeChecked(info.objectives.size(), true); + individualObjectivesToBeChecked.set(optimizingObjIndex, false); + do { //TODO introduce convergence criterion? (like max num of iterations) + WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); + refineSolution(separatingVector); + //Pick the threshold for the optimizing objective low enough so valid solutions are not excluded + for(auto const& paretoPoint : paretoOptimalPoints) { + thresholds[optimizingObjIndex] = std::min(thresholds[optimizingObjIndex], paretoPoint.first[optimizingObjIndex]); + } + if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ + std::cout << "INFEASIBLE OBJECTIVE THRESHOLDS" << std::endl; + return; + } + } while(!checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)); + STORM_LOG_DEBUG("Found a solution that satisfies the objective thresholds."); + individualObjectivesToBeChecked.clear(); + + // Improve the found solution. + // Note that we do not have to care whether a threshold is strict anymore, because the resulting optimum should be + // the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict + // thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds. + while (true) { + std::pair optimizationRes = underApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty."); + thresholds[optimizingObjIndex] = optimizationRes.first[optimizingObjIndex]; + STORM_LOG_DEBUG("Best solution found so far is " << storm::utility::convertNumber(thresholds[optimizingObjIndex]) << "."); + //Compute an upper bound for the optimum and check for convergence + optimizationRes = overApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + if(optimizationRes.second) { + RationalNumberType upperBoundOnOptimum = optimizationRes.first[optimizingObjIndex]; + STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber(upperBoundOnOptimum - thresholds[optimizingObjIndex])); + if(storm::utility::convertNumber(upperBoundOnOptimum - thresholds[optimizingObjIndex]) < 0.0001) { //TODO get this value from settings + std::cout << "Found Optimum: " << storm::utility::convertNumber(thresholds[optimizingObjIndex]) << "." << std::endl; + if(!checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)) { + std::cout << "Optimum is only the supremum" << std::endl; + } + return; + } + } + WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); + refineSolution(separatingVector); + } + } + + template + void SparseMultiObjectiveModelCheckerHelper::paretoQuery() { + //First consider the objectives individually + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + WeightVector direction(info.objectives.size(), storm::utility::zero()); + direction[objIndex] = storm::utility::one(); + refineSolution(direction); + } + + while(true) { //todo maxNumOfIterations ? + + // Get the halfspace of the underApproximation with maximal distance to a vertex of the overApproximation + std::vector> underApproxHalfspaces = underApproximation->getHalfspaces(); + std::vector overApproxVertices = overApproximation->getVertices(); + uint_fast64_t farestHalfspaceIndex = underApproxHalfspaces.size(); + RationalNumberType farestDistance = storm::utility::zero(); + for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < underApproxHalfspaces.size(); ++halfspaceIndex) { + for(auto const& vertex : overApproxVertices) { + RationalNumberType distance = underApproxHalfspaces[halfspaceIndex].euclideanDistance(vertex); + // Note that the distances are negative, i.e., we have to consider the minimum distance + if(distance < farestDistance) { + farestHalfspaceIndex = halfspaceIndex; + farestDistance = distance; + } + } + } + STORM_LOG_DEBUG("Farest distance between under- and overApproximation is " << storm::utility::convertNumber(farestDistance)); + if(storm::utility::convertNumber(farestDistance) > -0.0001) { //todo take this value from settings + std::cout << "DONE computing the pareto curve" << std::endl; + return; + } + WeightVector direction = underApproxHalfspaces[farestHalfspaceIndex].normalVector(); + // normalize the direction vector so the entries sum up to one + storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero())); + refineSolution(direction); + } + } + + template void SparseMultiObjectiveModelCheckerHelper::refineSolution(WeightVector const& direction) { weightedObjectivesChecker.check(storm::utility::vector::convertNumericVector(direction)); @@ -57,22 +192,7 @@ namespace storm { auto schedulerInsertRes = schedulers.insert(std::make_pair(std::move(scheduler), paretoOptimalPoints.end())); if(schedulerInsertRes.second){ // The scheduler is new, so insert the newly computed pareto optimal point. - // Due to numerical issues, however, it might be the case that the updated overapproximation will not contain the underapproximation, - // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case. - RationalNumberType computedOffset = storm::utility::vector::dotProduct(weightedObjectivesResult, direction); - RationalNumberType minimumOffset = computedOffset; - for(auto const& paretoPoint : paretoOptimalPoints){ - minimumOffset = std::max(minimumOffset, storm::utility::vector::dotProduct(paretoPoint.first, direction)); - } - if(minimumOffset > computedOffset){ - // We correct the issue by shifting the weightedObjectivesResult in the given direction such that - // it lies on the hyperplane given by direction and minimumOffset - WeightVector correction(direction); - storm::utility::vector::scaleVectorInPlace(correction, (minimumOffset - computedOffset) / storm::utility::vector::dotProduct(direction, direction)); - storm::utility::vector::addVectors(weightedObjectivesResult, correction, weightedObjectivesResult); - STORM_LOG_WARN("Numerical issues: The weightedObjectivesResult has to be shifted by " << storm::utility::sqrt(storm::utility::vector::dotProduct(correction, correction))); - } - // Now insert the point. For the corresponding scheduler, we safe the entry in the map. + // For the corresponding scheduler, we safe the entry in the map. schedulerInsertRes.first->second = paretoOptimalPoints.insert(std::make_pair(std::move(weightedObjectivesResult), std::vector())).first; } // In the case where the scheduler is not new, we assume that the corresponding pareto optimal points for the old and new scheduler are equal. @@ -88,17 +208,15 @@ namespace storm { } template - storm::storage::geometry::Halfspace SparseMultiObjectiveModelCheckerHelper::findSeparatingHalfspace(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) { + typename SparseMultiObjectiveModelCheckerHelper::WeightVector SparseMultiObjectiveModelCheckerHelper::findSeparatingVector(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) const { + STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by " << storm::utility::vector::convertNumericVector(pointToBeSeparated) << "."); // First, we check 'simple' weight vectors that correspond to checking a single objective while (!individualObjectivesToBeChecked.empty()){ uint_fast64_t objectiveIndex = individualObjectivesToBeChecked.getNextSetIndex(0); individualObjectivesToBeChecked.set(objectiveIndex, false); - WeightVector normalVector; // = (0..0 1 0..0) - normalVector.reserve(info.objectives.size()); - for(uint_fast64_t i = 0; i() : storm::utility::zero()); - } + WeightVector normalVector(info.objectives.size(), storm::utility::zero()); + normalVector[objectiveIndex] = storm::utility::one(); storm::storage::geometry::Halfspace h(normalVector, storm::utility::vector::dotProduct(normalVector, pointToBeSeparated)); bool hIsSeparating = true; @@ -106,10 +224,38 @@ namespace storm { hIsSeparating &= h.contains(paretoPoint.first); } if(hIsSeparating) { - return h; + return h.normalVector(); + } + } + + if(underApproximation->isEmpty()) { + // In this case [1 0..0] is always separating + WeightVector result(info.objectives.size(), storm::utility::zero()); + result.front() = storm::utility::one(); + return result; + } + + // Get the halfspace of the underApproximation with maximal distance to the given point + // Note that we are looking for a halfspace that does not contain the point. Thus the distances are negative and + // we are actually searching for the one with minimal distance. + std::vector> halfspaces = underApproximation->getHalfspaces(); + uint_fast64_t farestHalfspaceIndex = 0; + RationalNumberType farestDistance = halfspaces[farestHalfspaceIndex].euclideanDistance(pointToBeSeparated); + for(uint_fast64_t halfspaceIndex = 1; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) { + RationalNumberType distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); + if(distance < farestDistance) { + farestHalfspaceIndex = halfspaceIndex; + farestDistance = distance; } } + STORM_LOG_THROW(farestDistance<=storm::utility::zero(), storm::exceptions::UnexpectedException, "There is no seperating halfspace."); + // normalize the result so the entries sum up to one + WeightVector result = halfspaces[farestHalfspaceIndex].normalVector(); + storm::utility::vector::scaleVectorInPlace(result, storm::utility::one() / std::accumulate(result.begin(), result.end(), storm::utility::zero())); + return result; + + /* old version using LP solving // We now use LP solving to find a seperating halfspace. // Note that StoRM's LPSover does not support rational numbers. Hence, we optimize a polytope instead. // The variables of the LP are the normalVector entries w_0 ... w_{n-1} and the minimal distance d between the halfspace and the paretoOptimalPoints @@ -149,11 +295,24 @@ namespace storm { optimizeResult.first.pop_back(); //drop the distance RationalNumberType offset = storm::utility::vector::dotProduct(optimizeResult.first, pointToBeSeparated); return storm::storage::geometry::Halfspace(std::move(optimizeResult.first), std::move(offset)); + */ } template void SparseMultiObjectiveModelCheckerHelper::updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector) { storm::storage::geometry::Halfspace h(newWeightVector, storm::utility::vector::dotProduct(newWeightVector, newPoint)); + + // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation, + // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case. + RationalNumberType maximumOffset = h.offset(); + for(auto const& paretoPoint : paretoOptimalPoints){ + maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), paretoPoint.first)); + } + 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(newPoint)) << "."); + } overApproximation = overApproximation->intersection(h); STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); } @@ -167,8 +326,8 @@ namespace storm { } STORM_LOG_WARN("REMOVE ADDING ADDITIONAL VERTICES AS SOON AS HYPRO WORKS FOR DEGENERATED POLYTOPES"); - Point p1 = {-1, 0}; - Point p2 = {0, -1}; + Point p1 = {-1000, 0}; + Point p2 = {0, -1000}; paretoPointsVec.push_back(p1); paretoPointsVec.push_back(p2); @@ -191,6 +350,132 @@ namespace storm { STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); } + template + bool SparseMultiObjectiveModelCheckerHelper::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const { + std::vector> halfspaces = polytope->getHalfspaces(); + for(auto const& h : halfspaces) { + RationalNumberType 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; + } + } + } + } + return true; + } + + + + /* + + bool existsObjectiveWithoutThreshold = false; + bool existsObjectiveWithStrictThreshold = false; + Point thresholdPoint; + thresholdPoint.reserve(info.objectives.size()); + for(auto const& obj : info.objectives) { + if(obj.threshold) { + thresholdPoint.push_back(storm::utility::convertNumber(*obj.threshold)); + if(obj.thresholdIsStrict) { + existsObjectiveWithStrictThreshold = true; + } + } else { + existsObjectiveWithoutThreshold = true; + thresholdPoint.push_back(storm::utility::one()); + } + } + + if(existsObjectiveWithoutThreshold) { + // We need to find values for the objectives without thresholds in a way that the result is not affected by them. + // We build a polytope that contains valid solutions according to the thresholds and then analyze the intersection with the given polytope. + std::vector> thresholdConstraints; + WeightVector directionForOptimization(info.objectives.size(), storm::utility::zero()); + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + if(info.objectives[objIndex].threshold) { + WeightVector normalVector(info.objectives.size(), storm::utility::zero()); + normalVector[objIndex] = -storm::utility::one(); + thresholdConstraints.emplace_back(std::move(normalVector), -thresholdPoint[objIndex]); + directionForOptimization[objIndex] = -storm::utility::one(); + } + } + auto validSolutions = polytope->intersection(storm::storage::geometry::Polytope::create(std::move(thresholdConstraints))); + if(validSolutions->isEmpty()) { + return false; + } + auto optimizationRes = validSolutions->optimize(directionForOptimization); + STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "Couldn't check whether thresholds are satisfied for a given polytope: Didn't expect an infeasible or unbounded solution."); + + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + if(!info.objectives[objIndex].threshold) { + thresholdPoint[objIndex] = optimizationRes.first[objIndex] - storm::utility::one(); + } + } + } + + + + + + + // We then check whether the intersection of the thresholdPolytope and the given polytope contains any point (that also satisfies the strict thresholds); + + std::cout << "Valid solutions are " << validSolutions->toString(true) << std::endl; + std::cout << "not empty" << std::endl; + if(existsObjectiveWithStrictThreshold) { + std::cout << "hasObjWithStrictThreshold" << std::endl; + std::cout << "optRes is " << storm::utility::vector::convertNumericVector(optimizationRes.first) << std::endl; + for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + if(info.objectives[objIndex].threshold && info.objectives[objIndex].thresholdIsStrict) { + RationalNumberType threshold = storm::utility::convertNumber( *info.objectives[objIndex].threshold); + if (optimizationRes.first[objIndex] <= threshold) { + return false; + } + } + } + } + return true; + + + + Point thresholds; + thresholds.reserve(info.objectives.size()); + for(auto const& obj : info.objectives) { + if(obj.threshold) { + thresholds.push_back(storm::utility::convertNumber(*obj.threshold)); + } else { + STORM_LOG_ASSERT(numericalQueryResult, "Checking whether thresholds are satisfied for a numerical query but no threshold for the numerical objective was given."); + thresholds.push_back(*numericalQueryResult); + } + } + std::vector> halfspaces = polytope->getHalfspaces(); + for(auto const& h : halfspaces) { + RationalNumberType 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(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { + if(info.objectives[objIndex].thresholdIsStrict && h.normalVector()[objIndex] > storm::utility::zero()) { + return false; + } + } + } + } + return true; + + + + } + */ + #ifdef STORM_HAVE_CARL template class SparseMultiObjectiveModelCheckerHelper, storm::RationalNumber>; #endif diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h index 135eda82d..5de54b95a 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h @@ -24,6 +24,9 @@ namespace storm { typedef std::vector Point; typedef std::vector WeightVector; + static void check(Information& info); + + private: SparseMultiObjectiveModelCheckerHelper(Information& info); ~SparseMultiObjectiveModelCheckerHelper(); @@ -33,7 +36,6 @@ namespace storm { void paretoQuery(); - private: /* * Refines the solution w.r.t. the given direction vector @@ -41,17 +43,23 @@ namespace storm { void refineSolution(WeightVector const& direction); /* - * Returns a halfspace h that separates the underapproximation from the given point p, i.e., - * - p lies on the border of h and - * - for each x in the underApproximation, it holds that h contains x + * Returns a weight vector w that separates the underapproximation from the given point p, i.e., + * For each x in the underApproximation, it holds that w*x <= w*p * * @param pointToBeSeparated the point that is to be seperated * @param objectivesToBeCheckedIndividually stores for each objective whether it still makes sense to check for this objective individually (i.e., with weight vector given by w_{objIndex} = 1 ) */ - storm::storage::geometry::Halfspace findSeparatingHalfspace(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked); + WeightVector findSeparatingVector(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) const; + void updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector); void updateUnderApproximation(); + /* + * Returns true iff there is a point in the given polytope that satisfies the given thresholds. + * It is assumed that the given polytope contains the downward closure of its vertices. + */ + bool checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const; + Information& info; SparseWeightedObjectivesModelCheckerHelper weightedObjectivesChecker; diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h index c44f3c972..3b33227b1 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h @@ -23,6 +23,7 @@ namespace storm { std::shared_ptr originalFormula; bool originalFormulaMinimizes; boost::optional threshold; + bool thresholdIsStrict; std::string rewardModelName; boost::optional stepBound; @@ -32,9 +33,10 @@ namespace storm { out << (originalFormulaMinimizes ? "minimizes, \t" : "maximizes, \t"); out << "intern threshold:"; if(threshold){ + out << (thresholdIsStrict ? " >" : ">="); out << std::setw(5) << (*threshold) << ","; } else { - out << " none,"; + out << " none,"; } out << " \t"; out << "intern reward model: " << std::setw(10) << rewardModelName << ", \t"; diff --git a/src/storage/geometry/Halfspace.h b/src/storage/geometry/Halfspace.h index 0db210b0e..eb2bdbdcb 100644 --- a/src/storage/geometry/Halfspace.h +++ b/src/storage/geometry/Halfspace.h @@ -30,9 +30,29 @@ namespace storm { /* * Returns true iff the given point is contained in this halfspace, i.e., normalVector*point <= offset holds. */ - bool contains(std::vector const& point) { + bool contains(std::vector const& point) const { return storm::utility::vector::dotProduct(point, normalVector()) <= offset(); } + + /* + * 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 + * 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()); + } + + /* + * Returns the euclidean distance of the point from the hyperplane defined by this. + * The point is inside this halfspace iff 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 + return distance(point) / storm::utility::sqrt(storm::utility::vector::dotProduct(normalVector(), normalVector())); + } /* * Returns a string representation of this Halfspace.