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 <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::updateUnderApproximation(std::vector<RefinementStep> const& refinementSteps, std::shared_ptr<storm::storage::geometry::Polytope<RationalNumberType>>& underApproximation) {
- std::vector<Point> paretoPointsVec;
- paretoPointsVec.reserve(refinementSteps.size());
+ std::vector<Point> 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<Point> 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<RationalNumberType>::create(paretoPointsVec)->downwardClosure(upperBounds);
+ underApproximation = storm::storage::geometry::Polytope<RationalNumberType>::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<HyproPolytope<ValueType>>(internPolytope.linearTransformation(std::move(hyproMatrix), storm::adapters::toHypro(vector)));
}
- template <typename ValueType>
- std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::downwardClosure(boost::optional<Point> 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<HyproPolytope<ValueType>>(*this);
- }
- // Only keep the halfspaces where all entries of the normal vector are non-negative
- std::vector<hypro::Halfspace<ValueType>> halfspaces;
- for(auto& h : internPolytope.constraints()){
- if((h.normal().array() >= storm::utility::zero<ValueType>()).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<ValueType> direction = hypro::vector_t<ValueType>::Zero(internPolytope.dimension());
- direction(dim) = storm::utility::one<ValueType>();
- 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<ValueType> 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<HyproPolytope<ValueType>>(HyproPolytopeType(std::move(halfspaces)));
- }
-
-
template <typename ValueType>
std::pair<typename HyproPolytope<ValueType>::Point, bool> HyproPolytope<ValueType>::optimize(Point const& direction) const {
hypro::EvaluationResult<ValueType> 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<Polytope<ValueType>> linearTransformation(std::vector<Point> 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<Polytope<ValueType>> downwardClosure(boost::optional<Point> 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 <typename ValueType>
+ std::shared_ptr<Polytope<ValueType>> Polytope<ValueType>::createDownwardClosure(std::vector<Point> const& points) {
+ if(points.empty()) {
+ // In this case, the downwardclosure is empty
+ return createEmptyPolytope();
+ }
+ uint_fast64_t const dimensions = points.front().size();
+ std::vector<Halfspace<ValueType>> 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<Point> auxiliaryPoints = points;
+ auxiliaryPoints.reserve(auxiliaryPoints.size()*(1+dimensions));
+ for(auto const& point : points) {
+ for(uint_fast64_t dim=0; dim<dimensions; ++dim) {
+ auxiliaryPoints.push_back(point);
+ auxiliaryPoints.back()[dim] -= storm::utility::one<ValueType>();
+ }
+ }
+ std::vector<Halfspace<ValueType>> 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<ValueType>());
+ }
+ if(allGreaterZero){
+ halfspaces.push_back(std::move(h));
+ }
+ }
+ return create(halfspaces);
+ }
+
template <typename ValueType>
Polytope<ValueType>::Polytope() {
// Intentionally left empty
@@ -129,9 +166,8 @@ namespace storm {
}
template <typename ValueType>
- std::shared_ptr<Polytope<ValueType>> Polytope<ValueType>::downwardClosure(boost::optional<Point> const& upperBounds) const {
- STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented.");
- return nullptr;
+ std::shared_ptr<Polytope<ValueType>> Polytope<ValueType>::downwardClosure() const {
+ return createDownwardClosure(this->getVertices());
}
template <typename ValueType>
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<Polytope<ValueType>> 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<Polytope<ValueType>> createDownwardClosure(std::vector<Point> 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<Polytope<ValueType>> downwardClosure(boost::optional<Point> const& upperBounds = boost::none) const;
+ virtual std::shared_ptr<Polytope<ValueType>> downwardClosure() const;
/*!
* Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,