Browse Source

fix regarding creation of downward closure in 3D

Former-commit-id: f35406a83a
main
TimQu 9 years ago
parent
commit
669e9c6352
  1. 39
      src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
  2. 41
      src/storage/geometry/HyproPolytope.cpp
  3. 11
      src/storage/geometry/HyproPolytope.h
  4. 42
      src/storage/geometry/Polytope.cpp
  5. 14
      src/storage/geometry/Polytope.h

39
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));
}

41
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));

11
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).

42
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>

14
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.,

Loading…
Cancel
Save