TimQu
9 years ago
7 changed files with 547 additions and 69 deletions
-
15src/adapters/HyproAdapter.h
-
19src/storage/geometry/Halfspace.h
-
77src/storage/geometry/Hyperrectangle.h
-
195src/storage/geometry/HyproPolytope.cpp
-
124src/storage/geometry/HyproPolytope.h
-
106src/storage/geometry/Polytope.cpp
-
80src/storage/geometry/Polytope.h
@ -0,0 +1,77 @@ |
|||
#ifndef STORM_STORAGE_GEOMETRY_HYPERRECTANGLE_H_ |
|||
#define STORM_STORAGE_GEOMETRY_HYPERRECTANGLE_H_ |
|||
|
|||
#include <iostream> |
|||
#include <iomanip> |
|||
|
|||
#include "src/utility/macros.h" |
|||
#include "src/storage/geometry/Polytope.h" |
|||
#include "src/storage/geometry/Halfspace.h" |
|||
#include "src/exceptions/InvalidArgumentException" |
|||
|
|||
namespace storm { |
|||
namespace storage { |
|||
namespace geometry { |
|||
|
|||
/* |
|||
* This class represents a hyperrectangle, i.e., the intersection of finitely many intervals |
|||
*/ |
|||
|
|||
template <typename ValueType> |
|||
class Hyperrectangle { |
|||
|
|||
public: |
|||
|
|||
Hyperrectangle(std::vector<ValueType> const& lowerBounds, std::vector<ValueType> const& upperBounds) : mLowerBounds(lowerBounds), mUpperBounds(upperBounds) { |
|||
STORM_LOG_THROW(lowerBounds.size() == upperBounds.size(), storm::exceptions::InvalidArgumentException, "Tried to construct a hyperrectangle but the number of given lower bounds does not equal the number of given upper bounds."); |
|||
} |
|||
|
|||
Hyperrectangle(std::vector<ValueType>&& lowerBounds, std::vector<ValueType>&& upperBounds) : mLowerBounds(lowerBounds), mUpperBounds(upperBounds) { |
|||
STORM_LOG_THROW(lowerBounds.size() == upperBounds.size(), storm::exceptions::InvalidArgumentException, "Tried to construct a hyperrectangle but the number of given lower bounds does not equal the number of given upper bounds."); |
|||
} |
|||
|
|||
std::vector<ValueType> const& lowerBounds() const { |
|||
return mLowerBounds; |
|||
} |
|||
|
|||
std::vector<ValueType>& lowerBounds(){ |
|||
return mLowerBounds; |
|||
} |
|||
|
|||
std::vector<ValueType> const& upperBounds() const { |
|||
return mUpperBounds; |
|||
} |
|||
|
|||
std::vector<ValueType>& upperBounds(){ |
|||
return mUpperBounds; |
|||
} |
|||
|
|||
std::shared_ptr<Polytope<ValueType>> asPolytope() const { |
|||
STORM_LOG_THROW(lowerBounds.size() == upperBounds.size(), storm::exceptions::InvalidArgumentException, "Tried to construct a polytope form a hyperrectangle but the numbers of given lower and upper bounds do not match."); |
|||
std::vector<Halfspace<ValueType>> halfspaces; |
|||
halfspaces.reserve(2*lowerBounds.size()); |
|||
for(uint_fast64_t i = 0; i<lowerBounds.size(); ++i) { |
|||
std::vector<ValueType> direction(lowerBounds.size(), storm::utility::zero<ValueType>()); |
|||
direction[i] = -storm::utility::one<ValueType>(); |
|||
ValueType offset = -lowerBounds()[i]; |
|||
halfspaces.emplace_back(std::move(direction), std::move(offset)); |
|||
|
|||
direction = std::vector<ValueType>(lowerBounds.size(), storm::utility::zero<ValueType>()); |
|||
direction[i] = storm::utility::one<ValueType>(); |
|||
offset = upperBounds()[i]; |
|||
halfspaces.emplace_back(std::move(direction), std::move(offset)); |
|||
} |
|||
return Polytope<ValueType>::create(halfspaces); |
|||
} |
|||
|
|||
private: |
|||
|
|||
std::vector<ValueType> mLowerBounds; |
|||
std::vector<ValueType> mUpperBounds; |
|||
|
|||
}; |
|||
} |
|||
} |
|||
} |
|||
|
|||
#endif /* STORM_STORAGE_GEOMETRY_HYPERRECTANGLE_H_ */ |
@ -0,0 +1,195 @@ |
|||
#include "src/storage/geometry/HyproPolytope.h"
|
|||
|
|||
#include "src/utility/macros.h"
|
|||
#include "src/exceptions/InvalidArgumentException.h"
|
|||
#include "src/exceptions/UnexpectedException.h"
|
|||
|
|||
namespace storm { |
|||
namespace storage { |
|||
namespace geometry { |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(std::vector<Halfspace<ValueType>> const& halfspaces) { |
|||
if(halfspaces.empty()){ |
|||
internPolytope = HyproPolytopeType(); |
|||
} else { |
|||
std::vector<hypro::Halfspace<ValueType>> hyproHalfspaces; |
|||
hyproHalfspaces.reserve(halfspaces.size()); |
|||
for(auto& h : halfspaces) { |
|||
hyproHalfspaces.emplace_back(storm::adapters::toHypro(h)); |
|||
} |
|||
internPolytope = HyproPolytopeType(std::move(hyproHalfspaces)); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(std::vector<Point> const& points) { |
|||
if(points.empty()){ |
|||
internPolytope = HyproPolytopeType::Empty(); |
|||
} else { |
|||
std::vector<hypro::Point<ValueType>> hyproPoints; |
|||
hyproPoints.reserve(points.size()); |
|||
for(auto const& p : points){ |
|||
hyproPoints.emplace_back(storm::adapters::toHypro(p)); |
|||
} |
|||
internPolytope = HyproPolytopeType(std::move(hyproPoints)); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::create(boost::optional<std::vector<Halfspace<ValueType>>> const& halfspaces, |
|||
boost::optional<std::vector<Point>> const& points) { |
|||
if(halfspaces) { |
|||
STORM_LOG_WARN_COND(!points, "Creating a HyproPolytope where halfspaces AND points are given. The points will be ignored."); |
|||
return std::make_shared<HyproPolytope<ValueType>>(*halfspaces); |
|||
} else if(points) { |
|||
return std::make_shared<HyproPolytope<ValueType>>(*points); |
|||
} |
|||
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Creating a HyproPolytope but no representation was given."); |
|||
return nullptr; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(HyproPolytope<ValueType> const& other) : internPolytope(other.internPolytope) { |
|||
// Intentionally left empty
|
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(HyproPolytope<ValueType>&& other) : internPolytope(std::move(other.internPolytope)) { |
|||
// Intentionally left empty
|
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(HyproPolytopeType const& p) : internPolytope(p) { |
|||
// Intentionally left empty
|
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::HyproPolytope(HyproPolytopeType&& p) : internPolytope(p) { |
|||
// Intentionally left empty
|
|||
} |
|||
|
|||
template <typename ValueType> |
|||
HyproPolytope<ValueType>::~HyproPolytope() { |
|||
// Intentionally left empty
|
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::vector<typename Polytope<ValueType>::Point> HyproPolytope<ValueType>::getVertices() const { |
|||
std::vector<hypro::Point<ValueType>> hyproVertices = internPolytope.vertices(); |
|||
std::vector<Point> result; |
|||
result.reserve(hyproVertices.size()); |
|||
for(auto const& p : hyproVertices) { |
|||
result.push_back(storm::adapters::fromHypro(p.rawCoordinates())); |
|||
} |
|||
return result; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::vector<Halfspace<ValueType>> HyproPolytope<ValueType>::getHalfspaces() const { |
|||
std::vector<hypro::Halfspace<ValueType>> hyproHalfspaces = internPolytope.constraints(); |
|||
std::vector<Halfspace<ValueType>> result; |
|||
result.reserve(hyproHalfspaces.size()); |
|||
for(auto const& h : hyproHalfspaces) { |
|||
result.push_back(storm::adapters::fromHypro(h)); |
|||
} |
|||
return result; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
bool HyproPolytope<ValueType>::isEmpty() const { |
|||
return internPolytope.empty(); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
bool HyproPolytope<ValueType>::isUniversal() const { |
|||
return internPolytope.constraints().empty(); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
bool HyproPolytope<ValueType>::contains(Point const& point) const{ |
|||
return internPolytope.contains(storm::adapters::toHypro(point)); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
bool HyproPolytope<ValueType>::contains(std::shared_ptr<Polytope<ValueType>> const& other) const { |
|||
STORM_LOG_THROW(other->isHyproPolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a HyproPolytope and a different polytope implementation. This is not supported"); |
|||
return internPolytope.contains(other->asHyproPolytope().internPolytope); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::intersection(std::shared_ptr<Polytope<ValueType>> const& rhs) const{ |
|||
STORM_LOG_THROW(rhs->isHyproPolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a HyproPolytope and a different polytope implementation. This is not supported"); |
|||
return std::make_shared<HyproPolytope<ValueType>>(internPolytope.intersect(rhs->asHyproPolytope().internPolytope)); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::intersection(Halfspace<ValueType> const& halfspace) const{ |
|||
HyproPolytopeType result(internPolytope); |
|||
result.insert(storm::adapters::toHypro(halfspace)); |
|||
return std::make_shared<HyproPolytope<ValueType>>(std::move(result)); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::convexUnion(std::shared_ptr<Polytope<ValueType>> const& rhs) const{ |
|||
STORM_LOG_THROW(rhs->isHyproPolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a HyproPolytope and a different polytope implementation. This is not supported"); |
|||
return std::make_shared<HyproPolytope<ValueType>>(internPolytope.unite(rhs->asHyproPolytope().internPolytope)); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::minkowskiSum(std::shared_ptr<Polytope<ValueType>> const& rhs) const{ |
|||
STORM_LOG_THROW(rhs->isHyproPolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a HyproPolytope and a different polytope implementation. This is not supported"); |
|||
return std::make_shared<HyproPolytope<ValueType>>(internPolytope.minkowskiSum(rhs->asHyproPolytope().internPolytope)); |
|||
} |
|||
|
|||
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>>(internPolytope); |
|||
} |
|||
// 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(dim); |
|||
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> |
|||
bool HyproPolytope<ValueType>::isHyproPolytope() const { |
|||
return true; |
|||
} |
|||
|
|||
#ifdef STORM_HAVE_CARL
|
|||
template class HyproPolytope<storm::RationalNumber>; |
|||
#endif
|
|||
} |
|||
} |
|||
} |
@ -0,0 +1,124 @@ |
|||
#ifndef STORM_STORAGE_GEOMETRY_HYPROPOLYTOPE_H_ |
|||
#define STORM_STORAGE_GEOMETRY_HYPROPOLYTOPE_H_ |
|||
|
|||
#include "src/storage/geometry/Polytope.h" |
|||
#include "src/adapters/HyproAdapter.h" |
|||
|
|||
#ifdef STORM_HAVE_HYPRO |
|||
|
|||
namespace storm { |
|||
namespace storage { |
|||
namespace geometry { |
|||
|
|||
template <typename ValueType> |
|||
class HyproPolytope : public Polytope<ValueType> { |
|||
public: |
|||
|
|||
typedef typename Polytope<ValueType>::Point Point; |
|||
typedef hypro::HPolytope<ValueType> HyproPolytopeType; |
|||
|
|||
/*! |
|||
* Creates a HyproPolytope from the given halfspaces or points. |
|||
* If both representations are given, one of them might be ignored |
|||
*/ |
|||
static std::shared_ptr<Polytope<ValueType>> create(boost::optional<std::vector<Halfspace<ValueType>>> const& halfspaces, |
|||
boost::optional<std::vector<Point>> const& points); |
|||
|
|||
/*! |
|||
* Creates a HyproPolytope from the given halfspaces |
|||
* The resulting polytope is defined as the intersection of the halfspaces. |
|||
*/ |
|||
HyproPolytope(std::vector<Halfspace<ValueType>> const& halfspaces); |
|||
|
|||
/*! |
|||
* Creates a HyproPolytope from the given points. |
|||
* The resulting polytope is defined as the convex hull of the points' |
|||
*/ |
|||
HyproPolytope(std::vector<Point> const& points); |
|||
|
|||
/*! |
|||
* Copy and move constructors |
|||
*/ |
|||
HyproPolytope(HyproPolytope<ValueType> const& other); |
|||
HyproPolytope(HyproPolytope<ValueType>&& other); |
|||
|
|||
/*! |
|||
* Construction from a plain hypro polytope |
|||
*/ |
|||
HyproPolytope(HyproPolytopeType const& polytope); |
|||
HyproPolytope(HyproPolytopeType&& polytope); |
|||
|
|||
~HyproPolytope(); |
|||
|
|||
/*! |
|||
* Returns the vertices of this polytope. |
|||
*/ |
|||
virtual std::vector<Point> getVertices() const override; |
|||
|
|||
/*! |
|||
* Returns the halfspaces of this polytope. |
|||
*/ |
|||
virtual std::vector<Halfspace<ValueType>> getHalfspaces() const override; |
|||
|
|||
|
|||
/*! |
|||
* Returns whether this polytope is the empty set. |
|||
*/ |
|||
virtual bool isEmpty() const override; |
|||
|
|||
/*! |
|||
* Returns whether this polytope is universal (i.e., equals R^n). |
|||
*/ |
|||
virtual bool isUniversal() const override; |
|||
|
|||
/*! |
|||
* Returns true iff the given point is inside of the polytope. |
|||
*/ |
|||
virtual bool contains(Point const& point) const override; |
|||
|
|||
/*! |
|||
* Returns true iff the given polytope is a subset of this polytope. |
|||
*/ |
|||
virtual bool contains(std::shared_ptr<Polytope<ValueType>> const& other) const override; |
|||
|
|||
/*! |
|||
* Intersects this polytope with rhs and returns the result. |
|||
*/ |
|||
virtual std::shared_ptr<Polytope<ValueType>> intersection(std::shared_ptr<Polytope<ValueType>> const& rhs) const override; |
|||
virtual std::shared_ptr<Polytope<ValueType>> intersection(Halfspace<ValueType> const& halfspace) const override; |
|||
|
|||
/*! |
|||
* Returns the convex union of this polytope and rhs. |
|||
*/ |
|||
virtual std::shared_ptr<Polytope<ValueType>> convexUnion(std::shared_ptr<Polytope<ValueType>> const& rhs) const override; |
|||
|
|||
/*! |
|||
* Returns the minkowskiSum of this polytope and rhs. |
|||
*/ |
|||
virtual std::shared_ptr<Polytope<ValueType>> minkowskiSum(std::shared_ptr<Polytope<ValueType>> const& rhs) 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; |
|||
|
|||
virtual bool isHyproPolytope() const override; |
|||
|
|||
private: |
|||
|
|||
HyproPolytopeType internPolytope; |
|||
}; |
|||
|
|||
} |
|||
} |
|||
} |
|||
|
|||
#endif /* STORM_HAVE_HYPRO */ |
|||
|
|||
#endif /* STORM_STORAGE_GEOMETRY_HYPROPOLYTOPE_H_ */ |
Write
Preview
Loading…
Cancel
Save
Reference in new issue