From 669e9c6352dd2efae10279f72940d6e39e99056c Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 15 Jun 2016 18:42:18 +0200 Subject: [PATCH] fix regarding creation of downward closure in 3D Former-commit-id: f35406a83a88b903e526f2c31561d39a1eecbc8c --- .../helper/SparseMultiObjectiveHelper.cpp | 39 ++--------------- src/storage/geometry/HyproPolytope.cpp | 41 ------------------ src/storage/geometry/HyproPolytope.h | 11 ----- src/storage/geometry/Polytope.cpp | 42 +++++++++++++++++-- src/storage/geometry/Polytope.h | 14 +++---- 5 files changed, 50 insertions(+), 97 deletions(-) diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp index 62f48ce42..f2c04ea1d 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp @@ -243,43 +243,12 @@ namespace storm { template void SparseMultiObjectiveHelper::updateUnderApproximation(std::vector const& refinementSteps, std::shared_ptr>& underApproximation) { - std::vector paretoPointsVec; - paretoPointsVec.reserve(refinementSteps.size()); + std::vector paretoPoints; + paretoPoints.reserve(refinementSteps.size()); for(auto const& step : refinementSteps) { - paretoPointsVec.push_back(step.getPoint()); + paretoPoints.push_back(step.getPoint()); } - - STORM_LOG_WARN("REMOVE ADDING ADDITIONAL VERTICES AS SOON AS HYPRO WORKS FOR DEGENERATED POLYTOPES"); - if(paretoPointsVec.front().size()==2) { - Point p1 = {-10000, -9999}; - Point p2 = {-9999, -10000}; - paretoPointsVec.push_back(p1); - paretoPointsVec.push_back(p2); - } else { - Point p1 = {-10000, -9999, -9999}; - Point p2 = {-9999, -10000, -9999}; - Point p3 = {-9999, -9999, -10000}; - paretoPointsVec.push_back(p1); - paretoPointsVec.push_back(p2); - paretoPointsVec.push_back(p3); - } - - boost::optional upperBounds; - if(!paretoPointsVec.empty()){ - //Get the pointwise maximum of the pareto points - upperBounds = paretoPointsVec.front(); - for(auto paretoPointIt = paretoPointsVec.begin()+1; paretoPointIt != paretoPointsVec.end(); ++paretoPointIt){ - auto upperBoundIt = upperBounds->begin(); - for(auto const& paretoPointCoordinate : *paretoPointIt){ - if(paretoPointCoordinate>*upperBoundIt){ - *upperBoundIt = paretoPointCoordinate; - } - ++upperBoundIt; - } - } - } - - underApproximation = storm::storage::geometry::Polytope::create(paretoPointsVec)->downwardClosure(upperBounds); + underApproximation = storm::storage::geometry::Polytope::createDownwardClosure(paretoPoints); STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); } diff --git a/src/storage/geometry/HyproPolytope.cpp b/src/storage/geometry/HyproPolytope.cpp index 79c696a64..0670edd37 100644 --- a/src/storage/geometry/HyproPolytope.cpp +++ b/src/storage/geometry/HyproPolytope.cpp @@ -154,47 +154,6 @@ namespace storm { return std::make_shared>(internPolytope.linearTransformation(std::move(hyproMatrix), storm::adapters::toHypro(vector))); } - template - std::shared_ptr> HyproPolytope::downwardClosure(boost::optional const& upperBounds) const { - if(this->isUniversal() || this->isEmpty()) { - // In these cases, the polytope does not change, i.e., we return a copy of this - return std::make_shared>(*this); - } - // Only keep the halfspaces where all entries of the normal vector are non-negative - std::vector> halfspaces; - for(auto& h : internPolytope.constraints()){ - if((h.normal().array() >= storm::utility::zero()).all()){ - halfspaces.push_back(h); - } - } - // Add Halfspaces to bound the polytope in each (positive) direction - for(uint_fast64_t dim = 0; dim < internPolytope.dimension(); ++dim){ - hypro::vector_t direction = hypro::vector_t::Zero(internPolytope.dimension()); - direction(dim) = storm::utility::one(); - if(upperBounds){ - ValueType upperBound = (*upperBounds)[dim]; - halfspaces.emplace_back(std::move(direction), std::move(upperBound)); - } else { - // Compute max{x_dim | x \in P } - hypro::EvaluationResult evalRes = internPolytope.evaluate(direction); - switch (evalRes.errorCode) { - case hypro::SOLUTION::FEAS: - halfspaces.emplace_back(std::move(direction), std::move(evalRes.supportValue)); - break; - case hypro::SOLUTION::INFTY: - // if this polytope is not bounded in the current direction, no halfspace needs to be added - break; - default: - // there should never be an infeasible solution as we already excluded the empty polytope. - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected eror code for Polytope evaluation"); - break; - } - } - } - return std::make_shared>(HyproPolytopeType(std::move(halfspaces))); - } - - template std::pair::Point, bool> HyproPolytope::optimize(Point const& direction) const { hypro::EvaluationResult evalRes = internPolytope.evaluate(storm::adapters::toHypro(direction)); diff --git a/src/storage/geometry/HyproPolytope.h b/src/storage/geometry/HyproPolytope.h index 9de16fae8..78eb56ccd 100644 --- a/src/storage/geometry/HyproPolytope.h +++ b/src/storage/geometry/HyproPolytope.h @@ -106,17 +106,6 @@ namespace storm { */ virtual std::shared_ptr> linearTransformation(std::vector const& matrix, Point const& vector) const override; - /*! - * Returns the downward closure of this, i.e., the set { x | ex. y \in P : x<=y} where P is this Polytope. - * Put differently, the resulting polytope corresponds to this polytope, where - * 1. a vector y with y_i=max{x_i | x \in P} is computed and for each i, a halfspace with offset y_i and - * normal vector n (where n_i = 1 and the remaining entries are 0) is inserted. - * 2. all halfspaces where the normal vector has at least one negative entry are removed - * - * @param upperBounds If given, this vector is considered for y (hence, max{x_i | x i \in P does not need to be computed) - */ - virtual std::shared_ptr> downwardClosure(boost::optional const& upperBounds = boost::none) const override; - /*! * Finds an optimal point inside this polytope w.r.t. the given direction, i.e., * a point that maximizes dotPorduct(point, direction). diff --git a/src/storage/geometry/Polytope.cpp b/src/storage/geometry/Polytope.cpp index f5fc8ee6d..92f7b08c8 100644 --- a/src/storage/geometry/Polytope.cpp +++ b/src/storage/geometry/Polytope.cpp @@ -52,6 +52,43 @@ namespace storm { return nullptr; } + template + std::shared_ptr> Polytope::createDownwardClosure(std::vector const& points) { + if(points.empty()) { + // In this case, the downwardclosure is empty + return createEmptyPolytope(); + } + uint_fast64_t const dimensions = points.front().size(); + std::vector> halfspaces; + // We build the convex hull of the given points. + // However, auxiliary points (that will always be in the downward closure) are added. + // Then, the halfspaces of the resulting polytope are a superset of the halfspaces of the downward closure. + std::vector auxiliaryPoints = points; + auxiliaryPoints.reserve(auxiliaryPoints.size()*(1+dimensions)); + for(auto const& point : points) { + for(uint_fast64_t dim=0; dim(); + } + } + std::vector> auxiliaryHalfspaces = create(auxiliaryPoints)->getHalfspaces(); + // The downward closure is obtained by selecting the halfspaces for which the normal vector is non-negative (coordinate wise). + // Note that due to the auxiliary points, the polytope is never degenerated and thus there is always one unique halfspace-representation which is necessary: + // Consider, e.g., the convex hull of two points (1,0,0) and (0,1,1). + // There are multiple halfspace-representations for this set. In particular, there is one where all but one normalVectors have negative entries. + // However, the downward closure of this set can only be represented with 5 halfspaces. + for(auto& h : auxiliaryHalfspaces){ + bool allGreaterZero = true; + for(auto const& value : h.normalVector()) { + allGreaterZero &= (value >= storm::utility::zero()); + } + if(allGreaterZero){ + halfspaces.push_back(std::move(h)); + } + } + return create(halfspaces); + } + template Polytope::Polytope() { // Intentionally left empty @@ -129,9 +166,8 @@ namespace storm { } template - std::shared_ptr> Polytope::downwardClosure(boost::optional const& upperBounds) const { - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented."); - return nullptr; + std::shared_ptr> Polytope::downwardClosure() const { + return createDownwardClosure(this->getVertices()); } template diff --git a/src/storage/geometry/Polytope.h b/src/storage/geometry/Polytope.h index de08566e6..2950d2c94 100644 --- a/src/storage/geometry/Polytope.h +++ b/src/storage/geometry/Polytope.h @@ -41,6 +41,12 @@ namespace storm { */ static std::shared_ptr> createEmptyPolytope(); + /*! + * Creates the downward closure of the given points (i.e., the set { x | ex. y \in conv(points) : x<=y } + * If the vector of points is empty, the resulting polytope be empty. + */ + static std::shared_ptr> createDownwardClosure(std::vector const& points); + /*! * Returns the vertices of this polytope. */ @@ -98,14 +104,8 @@ namespace storm { /*! * Returns the downward closure of this, i.e., the set { x | ex. y \in P : x<=y} where P is this Polytope. - * Put differently, the resulting polytope corresponds to this polytope, where - * 1. a vector y with y_i=max{x_i | x \in P} is computed and for each i, a halfspace with offset y_i and - * normal vector n (where n_i = 1 and the remaining entries are 0) is inserted. - * 2. all halfspaces where the normal vector has at least one negative entry are removed - * - * @param upperBounds If given, this vector is considered for y (hence, max{x_i | x i \in P does not need to be computed) */ - virtual std::shared_ptr> downwardClosure(boost::optional const& upperBounds = boost::none) const; + virtual std::shared_ptr> downwardClosure() const; /*! * Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,