diff --git a/src/storm/storage/geometry/NativePolytope.cpp b/src/storm/storage/geometry/NativePolytope.cpp index 1f7c94567..2d35368d1 100644 --- a/src/storm/storage/geometry/NativePolytope.cpp +++ b/src/storm/storage/geometry/NativePolytope.cpp @@ -395,6 +395,11 @@ namespace storm { } return result; } + + template + std::shared_ptr> NativePolytope::clean() { + return create(boost::none, getVertices()); + } template class NativePolytope; diff --git a/src/storm/storage/geometry/NativePolytope.h b/src/storm/storage/geometry/NativePolytope.h index 11931a7d8..abf147bda 100644 --- a/src/storm/storage/geometry/NativePolytope.h +++ b/src/storm/storage/geometry/NativePolytope.h @@ -2,7 +2,6 @@ #define STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPE_H_ #include "storm/storage/geometry/Polytope.h" -#include "storm/storage/expressions/Expressions.h" #include "storm/adapters/EigenAdapter.h" namespace storm { @@ -125,8 +124,22 @@ namespace storm { */ virtual std::pair optimize(Point const& direction) const override; + /*! + * declares one variable for each dimension and returns the obtained variables. + * @param manager The expression manager that keeps track of the variables + * @param namePrefix The prefix that is prepanded to the variable index + */ + virtual std::vector declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const override; + + /*! + * returns the constrains defined by this polytope as an expression over the given variables + */ + virtual std::vector getConstraints(storm::expressions::ExpressionManager const& manager, std::vector const& variables) const override; + + virtual bool isNativePolytope() const override; + virtual std::shared_ptr> clean() override; private: // returns the vertices of this polytope as EigenVectors @@ -136,12 +149,6 @@ namespace storm { std::pair optimize(EigenVector const& direction) const; - // declares one variable for each constraint and returns the obtained variables. - std::vector declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const; - - // returns the constrains defined by this polytope as an expresseion - std::vector getConstraints(storm::expressions::ExpressionManager const& manager, std::vector const& variables) const; - //Stores whether the polytope is empty or not mutable EmptyStatus emptyStatus; diff --git a/src/storm/storage/geometry/Polytope.cpp b/src/storm/storage/geometry/Polytope.cpp index a89b7247e..1839e636f 100644 --- a/src/storm/storage/geometry/Polytope.cpp +++ b/src/storm/storage/geometry/Polytope.cpp @@ -47,35 +47,47 @@ namespace storm { template std::shared_ptr> Polytope::createDownwardClosure(std::vector const& points) { - if(points.empty()) { + if (points.empty()) { // In this case, the downwardclosure is empty return createEmptyPolytope(); } - uint_fast64_t const dimensions = points.front().size(); + // Reduce this call to a more general method + storm::storage::BitVector dimensions(points.front().size(), true); + return createSelectiveDownwardClosure(points, dimensions); + } + + template + std::shared_ptr> Polytope::createSelectiveDownwardClosure(std::vector const& points, storm::storage::BitVector const& selectedDimensions) { + if (points.empty()) { + // In this case, the downwardclosure is empty + return createEmptyPolytope(); + } + if (selectedDimensions.empty()) { + return create(points); + } + assert(points.front().size() == selectedDimensions.size()); + uint64_t dimensions = selectedDimensions.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. + // However, auxiliary points (that will always be in the selective 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()); + // The downward closure is obtained by erasing the halfspaces for which the normal vector is negative for one of the selected dimensions. + for (auto& h : auxiliaryHalfspaces) { + bool allGreaterEqZero = true; + for (auto const& dim : selectedDimensions) { + allGreaterEqZero &= (h.normalVector()[dim] >= storm::utility::zero()); } - if(allGreaterZero){ + if (allGreaterEqZero){ halfspaces.push_back(std::move(h)); } } @@ -155,11 +167,55 @@ namespace storm { return std::vector(); } + template + std::shared_ptr> Polytope::shift(Point const& b) const { + // perform an affine transformation with identity matrix + std::vector idMatrix(b.size(), Point(b.size(), storm::utility::zero())); + for (uint64_t i = 0; i < b.size(); ++i) { + idMatrix[i][i] = storm::utility::one(); + } + return affineTransformation(idMatrix, b); + } + + template + std::vector>> Polytope::setMinus(std::shared_ptr> const& rhs) const { + std::vector>> result; + auto rhsHalfspaces = rhs->getHalfspaces(); + std::shared_ptr> remaining = nullptr; + for (auto const& h : rhsHalfspaces) { + Polytope const& ref = (remaining == nullptr) ? *this : *remaining; + auto next = ref.intersection(h.invert()); + if (!next->isEmpty()) { + result.push_back(next); + } + remaining = ref.intersection(h); + if (remaining->isEmpty()) { + break; + } + } + return result; + } + template std::shared_ptr> Polytope::downwardClosure() const { return createDownwardClosure(this->getVertices()); } + template + std::vector Polytope::declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented for this polytope implementation."); + std::vector result; + return result; + } + + template + std::vector Polytope::getConstraints(storm::expressions::ExpressionManager const& manager, std::vector const& variables) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented for this polytope implementation."); + std::vector result; + return result; + } + + template template std::shared_ptr> Polytope::convertNumberRepresentation() const { @@ -197,6 +253,12 @@ namespace storm { return false; } + template + std::shared_ptr> Polytope::clean() { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "functionality not implemented for this polytope type."); + return nullptr; + } + template class Polytope; template std::shared_ptr> Polytope::convertNumberRepresentation() const; diff --git a/src/storm/storage/geometry/Polytope.h b/src/storm/storage/geometry/Polytope.h index d467af99d..590057215 100644 --- a/src/storm/storage/geometry/Polytope.h +++ b/src/storm/storage/geometry/Polytope.h @@ -113,6 +113,14 @@ namespace storm { */ virtual std::shared_ptr> downwardClosure() const; + /*! + * Computes the set {x \in this | x \notin rhs}. + * As this set is not necessarily convex, it is represented as the union of the returned polytopes. + * The returned polytopes are disjoint and non-empty. + * If the resulting set is empty, an empty vector is returned. + */ + std::vector>> setMinus(std::shared_ptr> const& rhs) const; + /*! * Finds an optimal point inside this polytope w.r.t. the given direction, i.e., * a point that maximizes dotPorduct(point, direction). @@ -128,6 +136,11 @@ namespace storm { template std::shared_ptr> convertNumberRepresentation() const; + /*! + * Performs cleaning operations, e.g., deleting redundant halfspaces + */ + virtual std::shared_ptr> clean(); + /* * Returns a string representation of this polytope. * @param numbersAsDouble If true, the occurring numbers are converted to double before printing to increase readability.