From 5cae7fca20fa312b08c50ea0bf21e5ec80596fe5 Mon Sep 17 00:00:00 2001 From: TimQu Date: Mon, 13 Feb 2017 15:27:51 +0100 Subject: [PATCH] started on native polytopes --- src/storm/adapters/EigenAdapter.cpp | 19 +- src/storm/adapters/EigenAdapter.h | 6 + src/storm/storage/geometry/NativePolytope.cpp | 290 ++++++++++ src/storm/storage/geometry/NativePolytope.h | 157 ++++++ .../HyperplaneCollector.cpp | 80 +++ .../HyperplaneCollector.h | 63 +++ .../HyperplaneEnumeration.cpp | 147 +++++ .../HyperplaneEnumeration.h | 68 +++ .../nativepolytopeconversion/QuickHull.cpp | 505 ++++++++++++++++++ .../nativepolytopeconversion/QuickHull.h | 165 ++++++ .../SubsetEnumerator.cpp | 93 ++++ .../SubsetEnumerator.h | 69 +++ 12 files changed, 1661 insertions(+), 1 deletion(-) create mode 100644 src/storm/storage/geometry/NativePolytope.cpp create mode 100644 src/storm/storage/geometry/NativePolytope.h create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp create mode 100644 src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h diff --git a/src/storm/adapters/EigenAdapter.cpp b/src/storm/adapters/EigenAdapter.cpp index 499488217..61cf7b98f 100644 --- a/src/storm/adapters/EigenAdapter.cpp +++ b/src/storm/adapters/EigenAdapter.cpp @@ -19,12 +19,29 @@ namespace storm { result->setFromTriplets(triplets.begin(), triplets.end()); return result; } - + + template + std::vector EigenAdapter::toStdVector(StormEigen::Matrix const& v){ + return std::vector(v.data(), v.data() + v.rows()); + } + + template + StormEigen::Matrix EigenAdapter::toEigenVector(std::vector const& v){ + return StormEigen::Matrix::Map(v.data(), v.size()); + } + template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + template std::vector EigenAdapter::toStdVector(StormEigen::Matrix const& v); + template StormEigen::Matrix EigenAdapter::toEigenVector(std::vector const& v); #ifdef STORM_HAVE_CARL template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + template std::vector EigenAdapter::toStdVector(StormEigen::Matrix const& v); + template StormEigen::Matrix EigenAdapter::toEigenVector(std::vector const& v); + template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + template std::vector EigenAdapter::toStdVector(StormEigen::Matrix const& v); + template StormEigen::Matrix EigenAdapter::toEigenVector(std::vector const& v); #endif } } diff --git a/src/storm/adapters/EigenAdapter.h b/src/storm/adapters/EigenAdapter.h index 6c5f8e7d5..a377afb87 100644 --- a/src/storm/adapters/EigenAdapter.h +++ b/src/storm/adapters/EigenAdapter.h @@ -17,6 +17,12 @@ namespace storm { */ template static std::unique_ptr> toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + + template + static std::vector toStdVector(StormEigen::Matrix const& v); + + template + static StormEigen::Matrix toEigenVector(std::vector const& v); }; } diff --git a/src/storm/storage/geometry/NativePolytope.cpp b/src/storm/storage/geometry/NativePolytope.cpp new file mode 100644 index 000000000..c45809551 --- /dev/null +++ b/src/storm/storage/geometry/NativePolytope.cpp @@ -0,0 +1,290 @@ +#include "storm/storage/geometry/NativePolytope.h" + +#ifdef STasdf + +#include "storm/utility/macros.h" +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/utility/solver.h" + +namespace storm { + namespace storage { + namespace geometry { + + template + NativePolytope::NativePolytope(std::vector> const& halfspaces) { + if(halfspaces.empty()){ + // The polytope is universal + emptyStatus = EmptyStatus::Nonempty; + } else { + uint_fast64_t maxCol = halfspaces.front().normalVector().size(); + uint_fast64_t maxRow = halfspaces.size(); + A = EigenMatrix(maxRow, maxCol); + b = EigenVector(maxRow); + for ( uint_fast64_t row = 0; row < maxRow; ++row ){ + assert(halfspaces[row].normal().rows() == maxCol); + b(row) = halfspaces[row].offset(); + A(row) = storm::adapters::EigenAdapter::toEigenVector(halfspaces[row].normal()); + } + emptyStatus = EmptyStatus::Unknown; + } + } + + template + NativePolytope::NativePolytope(std::vector const& points) { + if(points.empty()){ + emptyStatus = EmptyStatus::Empty; + } else { + std::vector eigenPoints; + eigenPoints.reserve(points.size()); + for(auto const& p : points){ + eigenPoints.emplace_back(storm::adapters::EigenAdapter::toEigenVector(p)); + } + + // todo: quickhull(eigenPoints) + + emptyStatus = EmptyStatus::Nonempty; + } + } + + template + std::shared_ptr> NativePolytope::create(boost::optional>> const& halfspaces, + boost::optional> const& points) { + if(halfspaces) { + STORM_LOG_WARN_COND(!points, "Creating a NativePolytope where halfspaces AND points are given. The points will be ignored."); + return std::make_shared>(*halfspaces); + } else if(points) { + return std::make_shared>(*points); + } + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Creating a NativePolytope but no representation was given."); + return nullptr; + } + + template + NativePolytope::NativePolytope(NativePolytope const& other) : emptyStatus(other.emptyStatus), A(other.A), b(other.b) { + // Intentionally left empty + } + + template + NativePolytope::NativePolytope(NativePolytope&& other) : emptyStatus(std::move(other.emptyStatus)), A(std::move(other.A)), b(std::move(other.b)) { + // Intentionally left empty + } + + template + NativePolytope::NativePolytope(EmptyStatus const& emptyStatus, EigenMatrix const& halfspaceMatrix, EigenVector const& halfspaceVector) : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) { + // Intentionally left empty + } + + template + NativePolytope::NativePolytope(EmptyStatus&& emptyStatus, EigenMatrix&& halfspaceMatrix, EigenVector&& halfspaceVector) : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) { + // Intentionally left empty + } + + template + NativePolytope::~NativePolytope() { + // Intentionally left empty + } + + template + std::vector::Point> NativePolytope::getVertices() const { + std::vector> eigenVertices = getEigenVertices(); + std::vector result; + result.reserve(eigenVertices.size()); + for(auto const& p : eigenVertices) { + result.push_back(storm::adapters::EigenAdapter::toStdVector(p)); + } + return result; + } + + template + std::vector> NativePolytope::getHalfspaces() const { + std::vector> result; + result.reserve(A.rows()); + + for(uint_fast64_t row=0; row < A.rows(); ++row){ + result.emplace_back(storm::adapters::EigenAdapter::toStdVector(A.row(row)), b.row(row)); + } + } + + template + bool NativePolytope::isEmpty() const { + if(emptyStatus == emptyStatus::Unknown) { + storm::expressions::ExpressionManager manager; + std::unique_ptr solver = storm::utility::solver::SmtSolverFactory().create(manager); + storm::expressions::Expression constraints = getConstraints(manager, declareVariables(manager, "x")); + solver->add(constraints); + switch(solver->check()) { + case storm::solver::SmtSolver::CheckResult::Sat: + emptyStatus = emptyStatus::Nonempty; + break; + case storm::solver::SmtSolver::CheckResult::Unsat: + emptyStatus = emptyStatus::Empty; + break; + default: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected result of SMT solver during emptyness-check of Polytope."); + break; + } + } + return emptyStatus == emptyStatus::Empty; + } + + template + bool NativePolytope::isUniversal() const { + return A.rows()==0; + } + + template + bool NativePolytope::contains(Point const& point) const{ + EigenVector x = storm::adapters::EigenAdapter::toEigenVector(point); + for(uint_fast64_t row=0; row < A.rows(); ++row){ + if((A.row(row) * x)(0) > b(row)){ + return false; + } + } + return true; + } + + template + bool NativePolytope::contains(std::shared_ptr> const& other) const { + STORM_LOG_THROW(other->isNativePolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported"); + // Check whether there is one point in other that is not in this + storm::expressions::ExpressionManager manager; + std::unique_ptr solver = storm::utility::solver::SmtSolverFactory().create(manager); + storm::expressions::Expression constraintsThis = this->getConstraints(manager, declareVariables(manager, "x")); + solver->add(!constraintsThis); + storm::expressions::Expression constraintsOther = dynamic_cast const&>(*other).getConstraints(manager, declareVariables(manager, "y")); + solver->add(constraintsOther); + + switch(solver->check()) { + case storm::solver::SmtSolver::CheckResult::Sat: + return false; + case storm::solver::SmtSolver::CheckResult::Unsat: + return true; + default: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected result of SMT solver during containment check of two polytopes."); + return false; + } + } + + template + std::shared_ptr> NativePolytope::intersection(std::shared_ptr> const& rhs) const{ + STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported"); + NativePolytope const& nativeRhs = dynamic_cast const&>(*rhs); + EigenMatrix resultA(A.rows() + nativeRhs.A.rows(), A.cols()); + resultA << A, + nativeRhs.A; + EigenVector resultb(resultA.rows()); + resultb << b, + nativeRhs.b; + return std::make_shared>(std::move(resultA), std::move(resultb)); + } + + template + std::shared_ptr> NativePolytope::intersection(Halfspace const& halfspace) const{ + EigenMatrix resultA(A.rows() + 1, A.cols()); + resultA << A, + storm::adapters::EigenAdapter::toEigenVector(halfspace.normal()); + EigenVector resultb(resultA.rows()); + resultb << b, + halfspace.offset(); + return std::make_shared>(std::move(resultA), std::move(resultb)); + } + + template + std::shared_ptr> NativePolytope::convexUnion(std::shared_ptr> const& rhs) const{ + STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported"); + if(this->isEmpty()) { + return std::make_shared>(dynamic_cast const&>(*rhs)); + } else if (rhs->isEmpty()) { + return std::make_shared>(*this); + } + if(this->isUniversal() || rhs->isUniversal) { + return std::make_shared>(*this); + } + + STORM_LOG_WARN("Implementation of convex union of two polytopes only works if the polytopes are bounded."); + + std::vector rhsVertices = dynamic_cast const&>(*rhs).getEigenVertices(); + std::vector resultVertices = this->getEigenVertices(); + resultVertices.insert(resultVertices.end(), std::make_move_iterator(rhsVertices.begin()), std::make_move_iterator(rhsVertices.end())); + + // todo invoke quickhull + + return nullptr; + } + + template + std::shared_ptr> NativePolytope::minkowskiSum(std::shared_ptr> const& rhs) const{ + STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException, "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported"); + NativePolytope const& nativeRhs = dynamic_cast const&>(*rhs); + + return std::make_shared>(internPolytope.minkowskiSum(dynamic_cast const&>(*rhs).internPolytope)); + } + + template + std::shared_ptr> NativePolytope::affineTransformation(std::vector const& matrix, Point const& vector) const{ + STORM_LOG_THROW(!matrix.empty(), storm::exceptions::InvalidArgumentException, "Invoked affine transformation with a matrix without rows."); + hypro::matrix_t hyproMatrix(matrix.size(), matrix.front().size()); + for(uint_fast64_t row = 0; row < matrix.size(); ++row) { + hyproMatrix.row(row) = storm::adapters::toNative(matrix[row]); + } + return std::make_shared>(internPolytope.affineTransformation(std::move(hyproMatrix), storm::adapters::toNative(vector))); + } + + template + std::pair::Point, bool> NativePolytope::optimize(Point const& direction) const { + hypro::EvaluationResult evalRes = internPolytope.evaluate(storm::adapters::toNative(direction)); + switch (evalRes.errorCode) { + case hypro::SOLUTION::FEAS: + return std::make_pair(storm::adapters::fromNative(evalRes.optimumValue), true); + case hypro::SOLUTION::UNKNOWN: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected eror code for Polytope evaluation"); + return std::make_pair(Point(), false); + default: + //solution is infinity or infeasible + return std::make_pair(Point(), false); + } + } + + template + bool NativePolytope::isNativePolytope() const { + return true; + } + template + std::vector::EigenVector> NativePolytope::getEigenVertices() const { + // todo: invoke conversion + return std::vector(); + } + + template + std::vector NativePolytope::declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const { + std::vector result; + result.reserve(A.cols()); + for(uint_fast64_t col=0; col < A.cols(); ++col){ + result.push_back(manager->declareVariable(namePrefix + std::to_string(col), manager->getRationalType())); + } + return result; + } + + template + storm::expressions::Expression NativePolytope::getConstraints(storm::expressions::ExpressionManager& manager, std::vector const& variables) cons { + storm::expressions::Expression result = manager.boolean(true); + for(uint_fast64_t row=0; row < A.rows(); ++row) { + storm::expressions::Expression lhs = manager.rational(A(row,0)) * variables[0].getExpression(); + for(uint_fast64_t col=1; col < A.cols(); ++col) { + lhs = lhs + manager.rational(A(row,col)) * variables[col].getExpression(); + } + result = result && (lhs <= manager.rational(b(row))); + } + return result; + } + + + template class NativePolytope; +#ifdef STORM_HAVE_CARL + template class NativePolytope; +#endif + } + } +} +#endif diff --git a/src/storm/storage/geometry/NativePolytope.h b/src/storm/storage/geometry/NativePolytope.h new file mode 100644 index 000000000..2155ac7ed --- /dev/null +++ b/src/storm/storage/geometry/NativePolytope.h @@ -0,0 +1,157 @@ +#ifndef STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPE_H_ +#define STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPE_H_ + +#include "storm/storage/geometry/Polytope.h" +#include "storm/storage/expressions/Expressions.h" +#include "storm/utility/eigen.h" + +#ifdef asdf + +namespace storm { + namespace storage { + namespace geometry { + + template + class NativePolytope : public Polytope { + public: + + typedef typename Polytope::Point Point; + + + /*! + * Creates a NativePolytope from the given halfspaces or points. + * If both representations are given, one of them might be ignored + */ + static std::shared_ptr> create(boost::optional>> const& halfspaces, + boost::optional> const& points); + + /*! + * Creates a NativePolytope from the given halfspaces + * The resulting polytope is defined as the intersection of the halfspaces. + */ + NativePolytope(std::vector> const& halfspaces); + + /*! + * Creates a NativePolytope from the given points. + * The resulting polytope is defined as the convex hull of the points' + */ + NativePolytope(std::vector const& points); + + /*! + * Copy and move constructors + */ + NativePolytope(NativePolytope const& other); + NativePolytope(NativePolytope&& other); + + /*! + * Construction from intern data + */ + NativePolytope(EmptyStatus const& emptyStatus, EigenMatrix const& halfspaceMatrix, EigenVector const& halfspaceVector); + NativePolytope(EmptyStatus&& emptyStatus, EigenMatrix&& halfspaceMatrix, EigenVector&& halfspaceVector); + + ~NativePolytope(); + + /*! + * Returns the vertices of this polytope. + */ + virtual std::vector getVertices() const override; + + /*! + * Returns the halfspaces of this polytope. + */ + virtual std::vector> 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> const& other) const override; + + /*! + * Intersects this polytope with rhs and returns the result. + */ + virtual std::shared_ptr> intersection(std::shared_ptr> const& rhs) const override; + virtual std::shared_ptr> intersection(Halfspace const& halfspace) const override; + + /*! + * Returns the convex union of this polytope and rhs. + */ + virtual std::shared_ptr> convexUnion(std::shared_ptr> const& rhs) const override; + + /*! + * Returns the minkowskiSum of this polytope and rhs. + */ + virtual std::shared_ptr> minkowskiSum(std::shared_ptr> const& rhs) const override; + + /*! + * Returns the affine transformation of this polytope P w.r.t. the given matrix A and vector b. + * The result is the set {A*x+b | x \in P} + * + * @param matrix the transformation matrix, given as vector of rows + * @param vector the transformation offset + */ + virtual std::shared_ptr> affineTransformation(std::vector const& matrix, Point const& vector) const override; + + /*! + * Finds an optimal point inside this polytope w.r.t. the given direction, i.e., + * a point that maximizes dotPorduct(point, direction). + * If such a point does not exist, the returned bool is false. There are two reasons for this: + * - The polytope is empty + * - The polytope is not bounded in the given direction + */ + virtual std::pair optimize(Point const& direction) const override; + + virtual bool isNativePolytope() const override; + + private: + typedef StormEigen::Matrix EigenMatrix; + typedef StormEigen::Matrix EigenVector; + + // returns the vertices of this polytope as EigenVectors + std::vector getEigenVertices() 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 + storm::expressions::Expression getConstraints(storm::expressions::ExpressionManager& manager, std::vector const& variables) const; + + + enum class EmptyStatus{ + Unknown, //It is unknown whether the polytope is empty or not + Empty, //The polytope is empty + Nonempty //the polytope is not empty + }; + //Stores whether the polytope is empty or not + mutable EmptyStatus emptyStatus; + + // Intern representation of the polytope as { x | Ax<=b } + EigenMatrix A; + EigenVector b; + + + + }; + + } + } +} + +#endif /* asdfsafO */ + +#endif /* STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPE_H_ */ diff --git a/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp new file mode 100644 index 000000000..395fe4d91 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp @@ -0,0 +1,80 @@ +/* + * File: HyperplaneCollector.cpp + * Author: tim quatmann + * + * Created on December 10, 2015, 6:06 PM + */ + +#include "HyperplaneCollector.h" +namespace hypro { + namespace pterm{ + + template< typename Number> + bool HyperplaneCollector::insert(hypro::Hyperplaneconst & hyperplane, std::vectorconst* indexList) { + return this->insert(hyperplane.normal(), hyperplane.offset(), indexList); + } + + template< typename Number> + bool HyperplaneCollector::insert(hypro::vector_t const& normal, Number const& offset, std::vectorconst* indexList) { + hypro::vector_t copyNormal(normal); + Number copyOffset(offset); + return this->insert(std::move(copyNormal), std::move(copyOffset), indexList); + } + + template< typename Number> + bool HyperplaneCollector::insert(hypro::vector_t && normal, Number && offset, std::vectorconst* indexList) { + //Normalize + Number infinityNorm = normal.template lpNorm(); + if(infinityNorm != (Number)0 ){ + normal /= infinityNorm; + offset /= infinityNorm; + } + + if(indexList == nullptr){ + //insert with empty list + return this->mMap.insert(ValueType(KeyType(normal, offset), std::vector())).second; + } else { + auto inserted = this->mMap.insert(ValueType(KeyType(normal, offset), *indexList)); + if(!inserted.second){ + //Append vertex list + inserted.first->second.insert(inserted.first->second.end(), indexList->begin(), indexList->end()); + } + return inserted.second; + } + } + + template< typename Number> + std::pair, hypro::vector_t> HyperplaneCollector::getCollectedHyperplanesAsMatrixVector() const{ + assert(!this->mMap.empty()); + hypro::matrix_t A(this->mMap.size(), this->mMap.begin()->first.first.rows()); + hypro::vector_t b(this->mMap.size()); + + std::size_t row = 0; + for(auto const& mapEntry : this->mMap){ + A.row(row) = mapEntry.first.first; + b(row) = mapEntry.first.second; + ++row; + } + return std::pair, hypro::vector_t>(std::move(A), std::move(b)); + } + + template< typename Number> + std::vector> HyperplaneCollector::getIndexLists() const{ + assert(!this->mMap.empty()); + std::vector> result(this->mMap.size()); + + auto resultIt = result.begin(); + for(auto const& mapEntry : this->mMap){ + *resultIt = mapEntry.second; + ++resultIt; + } + return result; + } + + template< typename Number> + std::size_t HyperplaneCollector::numOfCollectedHyperplanes() const { + return this->mMap.size(); + } + } +} + diff --git a/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h new file mode 100644 index 000000000..d4a3c7d6f --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h @@ -0,0 +1,63 @@ +/* + * File: HyperplaneCollector.h + * Author: tim quatmann + * + * Created on December 10, 2015, 6:06 PM + */ + +#pragma once +#include +#include "../../config.h" +#include "../../datastructures/Hyperplane.h" + +namespace hypro{ + namespace pterm{ + /*! + * This class can be used to collect a set of hyperplanes (without duplicates). + * The inserted hyperplanes are normalized, i.e. devided by the infinity norm of the normal vector + */ + template< typename Number> + class HyperplaneCollector { + public: + HyperplaneCollector() = default; + HyperplaneCollector(const HyperplaneCollector& orig) = default; + virtual ~HyperplaneCollector() = default; + + /* + * inserts the given hyperplane. + * For every (unique) hyperplane, there is a list of indices which can be used e.g. to obtain the set of vertices that lie on each hyperplane. + * If indexList is given (i.e. not nullptr), the given indices are appended to that list. + * Returns true iff the hyperplane was inserted (i.e. the hyperplane was not already contained in this) + */ + bool insert(hypro::Hyperplane const& hyperplane, std::vectorconst* indexList = nullptr); + bool insert(hypro::vector_t const& normal, Number const& offset, std::vectorconst* indexList = nullptr); + bool insert(hypro::vector_t&& normal, Number && offset, std::vectorconst* indexList = nullptr); + + std::pair, hypro::vector_t> getCollectedHyperplanesAsMatrixVector() const; + //Note that the returned lists might contain dublicates. + std::vector> getIndexLists() const; + + std::size_t numOfCollectedHyperplanes() const; + + private: + typedef std::pair, Number> NormalOffset; + class NormalOffsetHash{ + public: + std::size_t operator()(NormalOffset const& ns) const { + std::size_t seed = std::hash>()(ns.first); + carl::hash_add(seed, std::hash()(ns.second)); + return seed; + } + }; + typedef typename std::unordered_map, NormalOffsetHash>::key_type KeyType; + typedef typename std::unordered_map, NormalOffsetHash>::value_type ValueType; + + std::unordered_map, NormalOffsetHash> mMap; + + + }; + } +} + + +#include "HyperplaneCollector.tpp" \ No newline at end of file diff --git a/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp new file mode 100644 index 000000000..5339afa24 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp @@ -0,0 +1,147 @@ +/* + * File: HyperplaneEnumeration.tpp + * Author: tim quatmann + * Author: phillip florian + * + * Created on December 28, 2015, 1:06 PM + */ + +#include "HyperplaneEnumeration.h" + +#include "../SubsetEnumerator.h" +#include "../HyperplaneCollector.h" + +namespace hypro { + namespace pterm{ + + template + bool HyperplaneEnumeration::generateVerticesFromHalfspaces(PTermHPolytope const& hPoly, bool generateRelevantHyperplanesAndVertexSets){ + PTERM_DEBUG("Invoked generateVerticesFromHalfspaces with " << hPoly.getMatrix().rows() << " hyperplanes. Dimension is " << hPoly.dimension()); + std::unordered_map, std::set> vertexCollector; + hypro::pterm::SubsetEnumerator> subsetEnum(hPoly.getMatrix().rows(), hPoly.dimension(), hPoly.getMatrix(), this->linearDependenciesFilter); + if(subsetEnum.setToFirstSubset()){ + do{ + std::vector const& subset = subsetEnum.getCurrentSubset(); + std::pair, hypro::vector_t> subHPolytope(hPoly.getSubHPolytope(subset)); + Point point(subHPolytope.first.fullPivLu().solve(subHPolytope.second)); + //Point point(hypro::gauss(subHPolytope.first, subHPolytope.second)); + if(hPoly.contains(point)){ + //Note that the map avoids duplicates. + auto hyperplaneIndices = vertexCollector.insert(typename std::unordered_map, std::set>::value_type(std::move(point), std::set())).first; + if(generateRelevantHyperplanesAndVertexSets){ + hyperplaneIndices->second.insert(subset.begin(), subset.end()); + } + } + } while (subsetEnum.incrementSubset()); + } else { + std::cout << "Could not generate any vertex while converting from H Polytope. TODO: implement this case (we get some unbounded thing here)" << std::endl; + return false; + } + if(generateRelevantHyperplanesAndVertexSets){ + //For each hyperplane, get the number of (unique) vertices that lie on it. + std::vector verticesOnHyperplaneCounter(hPoly.getMatrix().rows(), 0); + for(auto const& mapEntry : vertexCollector){ + for(auto const& hyperplaneIndex : mapEntry.second){ + ++verticesOnHyperplaneCounter[hyperplaneIndex]; + } + } + + //Only keep the hyperplanes on which at least dimension() vertices lie. + //Note that this will change the indices of the Hyperplanes. + //Therefore, we additionally store the old indices for every hyperplane to be able to translate from old to new indices + hypro::pterm::HyperplaneCollector hyperplaneCollector; + for(std::size_t hyperplaneIndex = 0; hyperplaneIndex < verticesOnHyperplaneCounter.size(); ++hyperplaneIndex){ + if(verticesOnHyperplaneCounter[hyperplaneIndex] >= hPoly.dimension()){ + std::vector oldIndex; + oldIndex.push_back(hyperplaneIndex); + hyperplaneCollector.insert(hPoly.getMatrix().row(hyperplaneIndex), hPoly.getVector()(hyperplaneIndex), &oldIndex); + } + } + auto matrixVector = hyperplaneCollector.getCollectedHyperplanesAsMatrixVector(); + this->mRelevantMatrix = std::move(matrixVector.first); + this->mRelevantVector = std::move(matrixVector.second); + + //Get the mapping from old to new indices + std::vector oldToNewIndexMapping (hPoly.getMatrix().rows(), hPoly.getMatrix().rows()); //Initialize with some illegal value + std::vector> newToOldIndexMapping(hyperplaneCollector.getIndexLists()); + for(std::size_t newIndex = 0; newIndex < newToOldIndexMapping.size(); ++newIndex){ + for(auto const& oldIndex : newToOldIndexMapping[newIndex]){ + oldToNewIndexMapping[oldIndex] = newIndex; + } + } + + //Insert the resulting vertices and get the set of vertices that lie on each hyperplane + std::vector> vertexSets(this->mRelevantMatrix.rows()); + this->mResultVertices.clear(); + this->mResultVertices.reserve(vertexCollector.size()); + for(auto const& mapEntry : vertexCollector){ + for(auto const& oldHyperplaneIndex : mapEntry.second){ + //ignore the hyperplanes which are redundant, i.e. for which there is no new index + if(oldToNewIndexMapping[oldHyperplaneIndex] < this->mRelevantVector.rows()){ + vertexSets[oldToNewIndexMapping[oldHyperplaneIndex]].insert(this->mResultVertices.size()); + } + } + this->mResultVertices.push_back(mapEntry.first); + } + this->mVertexSets.clear(); + this->mVertexSets.reserve(vertexSets.size()); + for(auto const& vertexSet : vertexSets){ + this->mVertexSets.emplace_back(vertexSet.begin(), vertexSet.end()); + } + + } else { + this->mResultVertices.clear(); + this->mResultVertices.reserve(vertexCollector.size()); + for(auto const& mapEntry : vertexCollector){ + this->mResultVertices.push_back(mapEntry.first); + } + this->mVertexSets.clear(); + this->mRelevantMatrix = hypro::matrix_t(); + this->mRelevantVector = hypro::vector_t(); + } + PTERM_DEBUG("Invoked generateVerticesFromHalfspaces with " << hPoly.getMatrix().rows() << " hyperplanes and " << this->mResultVertices.size() << " vertices and " << this->mRelevantMatrix.rows() << " relevant hyperplanes. Dimension is " << hPoly.dimension()); + return true; + } + + + template + bool HyperplaneEnumeration::linearDependenciesFilter(std::vector const& subset, std::size_t const& item, hypro::matrix_t const& A) { + hypro::matrix_t subMatrix(subset.size() +1, A.cols()); + for (std::size_t i = 0; i < subset.size(); ++i){ + subMatrix.row(i) = A.row(subset[i]); + } + subMatrix.row(subset.size()) = A.row(item); + Eigen::FullPivLU> lUMatrix( subMatrix ); + // std::cout << "The rank is " << lUMatrix.rank() << " and the matrix has " << subMatrix.rows() << " rows" << std::endl; + if (lUMatrix.rank() < subMatrix.rows()){ + //Linear dependent! + return false; + } else { + return true; + } + } + + + template + std::vector>& HyperplaneEnumeration::getResultVertices() { + return this->mResultVertices; + } + + template + hypro::matrix_t& HyperplaneEnumeration::getRelevantMatrix() { + return this->mRelevantMatrix; + } + + template + hypro::vector_t& HyperplaneEnumeration::getRelevantVector() { + return this->mRelevantVector; + } + + template + std::vector>& HyperplaneEnumeration::getVertexSets() { + return this->mVertexSets; + } + + } +} + diff --git a/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h new file mode 100644 index 000000000..3ba77c1d2 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h @@ -0,0 +1,68 @@ +/* + * File: HyperplaneEnumeration.h + * Author: tim quatmann + * Author: phillip florian + * + * Created on December 27, 2015, 1:06 PM + */ + +#pragma once +#include "../macros.h" +#include "../../../datastructures/Hyperplane.h" +#include "../../../datastructures/Point.h" + +namespace hypro{ + namespace pterm{ + + template< typename Number> + class HyperplaneEnumeration { + public: + HyperplaneEnumeration() = default; + virtual ~HyperplaneEnumeration() = default; + + /* + * Generates the vertices of the given polytope by enumerating all intersection points generated by subsets of hyperplanes of size hPoly.dimension(). + * If the given flag is true, this method will also compute + * * the minimal set of hyperplanes which represent the given hPoly (can be used to remove redundant hyperplanes), and + * * for each hyperplane, the set of (non-redundant) vertices that lie on that hyperplane. + * + * Use the provided getter methods to retrieve the results + * + * @return true iff conversion was successful. + */ + bool generateVerticesFromHalfspaces(PTermHPolytope const& hPoly, bool generateRelevantHyperplanesAndVertexSets); + + std::vector>& getResultVertices(); + + /*! + * Returns the set of halfspaces which are not redundant + * @note the returned matrix and vector are empty if the corresponding flag was false + */ + hypro::matrix_t& getRelevantMatrix(); + hypro::vector_t& getRelevantVector(); + + /*! + * Returns for each hyperplane the set of vertices that lie on that hyperplane. + * A vertex is given as an index in the relevantVertices vector. + * @note the returned vector is empty if the corresponding flag was false + */ + std::vector>& getVertexSets(); + + + /* + * Returns true if the hyperplanes with indices of subset and item are all linear independent + * Note that this is also used by the hybrid polytope. + */ + static bool linearDependenciesFilter(std::vector const& subset, std::size_t const& item, hypro::matrix_t const& A); + + private: + std::vector> mResultVertices; + hypro::matrix_t mRelevantMatrix; + hypro::vector_t mRelevantVector; + std::vector> mVertexSets; + }; + } +} + + +#include "HyperplaneEnumeration.tpp" \ No newline at end of file diff --git a/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp b/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp new file mode 100644 index 000000000..f252925e9 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp @@ -0,0 +1,505 @@ +#include "storm/storage/geometry/nativepolytopeconversion/QuickHull.h" + +#include + +#include "storm/utility/macros.h" + +#include "storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h" +#include "storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h" + +namespace hypro { + namespace pterm{ + + template + void QuickHull::generateHalfspacesFromVertices(std::vector const& points, bool generateRelevantVerticesAndVertexSets){ + STORM_LOG_ASSERT(!points.empty(), "Invoked QuickHull with empty set of points."); + STORM_LOG_DEBUG("Invoked QuickHull on " << points.size() << " points"); + const uint_fast64_t dimension = points.front().rows(); + + // Generate initial set of d+1 affine independent points (if such a set exists) + std::vector vertexIndices; + uint_fast64_t minMaxVertexNumber; + if(!this->findInitialVertices(points, vertexIndices, minMaxVertexNumber)) { + // todo deal with this + std::cout << "QuickHull: Could not find d+1 affine independend points. TODO: implement this case (we get some degenerated thing here)" << std::endl; + } + + // compute point inside initial facet + EigenVector insidePoint(EigenVector::Zero(dimension)); + for(uint_fast64_t vertexIndex : vertexIndices){ + insidePoint += points[vertexIndex]; + } + insidePoint /= storm::utility::convertNumber(vertexIndices.size()); + + // Create the initial facets from the found vertices. + std::vector facets = computeInitialFacets(vertexIndices, insidePoint); + + // Enlarge the mesh by adding by first considering all points that are min (or max) in at least one dimension + storm::storage::BitVector currentFacets(facets.size(), true); + storm::storage::BitVector consideredPoints(points.size(), false); + uint_fast64_t currentNumOfVertices = vertexIndices.size(); + for(uint_fast64_t i = 0; i < minMaxVertexNumber; ++i) { + consideredPoints.set(i); + } + this->extendMesh(points, consideredPoints, facets, currentFacets, insidePoint, currentNumOfVertices); + for(auto & facet : facets){ + facet.maxOutsidePointIndex = 0; + facet.outsideSet.clear(); + } + + consideredPoints = storm::storage::BitVector(points.size(), true); + this->extendMesh(points, consideredPoints, facets, currentFacets, insidePoint, currentNumOfVertices); + this->getPolytopeFromMesh(points, facets, currentFacets, success && generateRelevantVerticesAndVertexSets); + + return true; + } + + + template + void QuickHull::extendMesh(std::vector& points, + storm::storage::BitVector& consideredPoints, + std::vector& facets, + storm::storage::BitVector& currentFacets, + vector_t& insidePoint, + uint_fast64_t& currentNumOfVertices) const { + + storm::storage::BitVector currentOutsidePoints = consideredPoints; + // Compute initial outside Sets + for(uint_fast64_t facetIndex : currentFacets){ + computeOutsideSetOfFacet(facets[facetIndex], currentOutsidePoints, points); + } + + for(uint_fast64_t facetCount = currentFacets.getNextSetIndex(0); facetCount != currentFacets.size(); facetCount = currentFacets.getNextSetIndex(facetCount+1)) { + // set all points to false to get rid of points that lie within the polytope after each iteration + currentOutsidePoints.reset(); + // Find a facet with a non-empty outside set + if(!facets[facetCount].outsideSet.empty()) { + uint_fast64_t numberOfNewFacets = 0; + // Now we compute the enlarged mesh + uint_fast64_t farAwayPointIndex = facets[facetCount].maxOutsidePointIndex; + assert(consideredPoints.get(farAwayPointIndex)); + // get Visible set from maxOutsidePoint of the current facet + std::set visibleSet = getVisibleSet(facets, facetCount, points[farAwayPointIndex]); + std::set invisibleSet = getInvisibleNeighbors(facets, visibleSet); + for(auto invisFacetIt = invisibleSet.begin(); invisFacetIt != invisibleSet.end(); ++invisFacetIt) { + for(auto visFacetIt = visibleSet.begin(); visFacetIt != visibleSet.end(); ++visFacetIt) { + if (facetHasNeighborWithIndex(facets[*invisFacetIt], *visFacetIt)) { + Facet newFacet; + // Set points of Facet + newFacet.points = getCommonPoints(facets[*invisFacetIt], facets[*visFacetIt]); + // replace old facet index by new facet index in the current neighbor + newFacet.points.push_back(farAwayPointIndex); + replaceFacetNeighbor(facets, *visFacetIt, facets.size(), *invisFacetIt); + newFacet.neighbors.push_back(*invisFacetIt); + // get new hyperplane from common points and new point + std::vector vectorSet; + vectorSet.reserve(points[0].dimension()); + EigenVector refPoint(points[farAwayPointIndex].rawCoordinates()); + for (uint_fast64_t pointCount = 0; pointCount + 1 < points[0].dimension(); ++pointCount){ + vectorSet.emplace_back(points[newFacet.points[pointCount]].rawCoordinates() - refPoint); + } + assert(linearDependenciesFilter(vectorSet)); + newFacet.hyperplane = Hyperplane(std::move(refPoint), std::move(vectorSet)); + // check orientation of hyperplane + // To avoid multiple matrix multiplications we need a point that is known to lie well within the polytope + if (!newFacet.hyperplane.holds(insidePoint)){ + newFacet.hyperplane.invert(); + } + facets.push_back(newFacet); + // avoid using replaced facets and add new ones + currentFacets.push_back(true); + // Set Points in outsideSet free + // increase Number Of new Facets + ++numberOfNewFacets; + } + } + } + + for(auto visFacetIt = visibleSet.begin(); visFacetIt != visibleSet.end(); ++visFacetIt){ + for(uint_fast64_t m = 0; m < facets[*visFacetIt].outsideSet.size(); ++m){ + currentOutsidePoints.set(facets[*visFacetIt].outsideSet[m], true); + } + } + for(auto visFacetIt = visibleSet.begin(); visFacetIt != visibleSet.end(); ++visFacetIt){ + currentFacets.set(*visFacetIt, false); + } + // compute new outside sets + for(uint_fast64_t fIndex = facets.size()-numberOfNewFacets; fIndex != facets.size(); ++fIndex){ + computeOutsideSetOfFacet(facets[fIndex], currentOutsidePoints, points); + } + + ++currentNumOfVertices; + +#ifdef PTERM_DEBUG_OUTPUT + numOfIncludedPoints += currentOutsidePoints.count(); + PTERM_DEBUG("Mesh currently contains " << numOfIncludedPoints << " of " << consideredPoints.count() << "points"); +#endif + + // find neighbors in new facets + setNeighborhoodOfNewFacets(facets, facets.size() - numberOfNewFacets, points[0].dimension()); + + } + } + return true; + } + + template + void QuickHull::getPolytopeFromMesh(std::vector const& points, std::vector const& facets, storm::storage::BitVector const& currentFacets, bool generateRelevantVerticesAndVertexSets){ + + hypro::pterm::HyperplaneCollector hyperplaneCollector; + for(uint_fast64_t facetCount = 0; facetCount < facets.size(); ++facetCount){ + if(currentFacets[facetCount]){ + hyperplaneCollector.insert(std::move(facets[facetCount].hyperplane), generateRelevantVerticesAndVertexSets ? &facets[facetCount].points : nullptr); + } + } + + if(generateRelevantVerticesAndVertexSets){ + //Get the mapping from a hyperplane to the set of vertices that lie on that plane, erase the duplicates, and count for each vertex the number of hyperplanes on which that vertex lies + this->mVertexSets = hyperplaneCollector.getIndexLists(); + std::vector hyperplanesOnVertexCounter(points.size(), 0); + for(auto& vertexVector : this->mVertexSets){ + std::set vertexSet; + for(auto const& i : vertexVector){ + if(vertexSet.insert(i).second){ + ++hyperplanesOnVertexCounter[i]; + } + } + vertexVector.assign( vertexSet.begin(), vertexSet.end()); + } + //Now, we can erase all vertices which do not lie on at least dimension() hyperplanes. + //Note that the indices of the HyperplaneToVerticesMapping needs to be refreshed according to the new set of vertices + //Therefore, we additionally store the old indices for every vertex to be able to translate from old to new indices + std::unordered_map> relevantVerticesMap; + relevantVerticesMap.reserve(points.size()); + for(uint_fast64_t vertexIndex = 0; vertexIndex < hyperplanesOnVertexCounter.size(); ++vertexIndex){ + if(hyperplanesOnVertexCounter[vertexIndex] >= points[0].dimension()){ + auto mapEntry = relevantVerticesMap.insert(typename std::unordered_map>::value_type(points[vertexIndex], std::vector())).first; + mapEntry->second.push_back(vertexIndex); + } + } + //Fill in the relevant vertices and create a translation map from old to new indices + std::vector oldToNewIndexMapping (points.size(), points.size()); //Initialize with some illegal value + this->mRelevantVertices.clear(); + this->mRelevantVertices.reserve(relevantVerticesMap.size()); + for(auto const& mapEntry : relevantVerticesMap){ + for(auto const& oldIndex : mapEntry.second){ + oldToNewIndexMapping[oldIndex] = this->mRelevantVertices.size(); + } + this->mRelevantVertices.push_back(mapEntry.first); + } + //Actually translate and erase duplicates + for(auto& vertexVector : this->mVertexSets){ + std::set vertexSet; + for(auto const& oldIndex : vertexVector){ + if(hyperplanesOnVertexCounter[oldIndex] >= points[0].dimension()){ + vertexSet.insert(oldToNewIndexMapping[oldIndex]); + } + } + vertexVector.assign( vertexSet.begin(), vertexSet.end()); + } + } else { + this->mRelevantVertices.clear(); + this->mVertexSets.clear(); + } + auto matrixVector = hyperplaneCollector.getCollectedHyperplanesAsMatrixVector(); + this->mResultMatrix = std::move(matrixVector.first); + this->mResultVector = std::move(matrixVector.second); + + PTERM_DEBUG("Computed H representation from " << points.size() << " vertices and " << this->mResultVector.rows() << " hyperplanes and " << this->mRelevantVertices.size() << " relevant vertices. Dimension is " << points[0].dimension()); + PTERM_DEBUG("Total number of considered facets: " << facets.size() << " where " << currentFacets.count() << " are enabled."); + } + + template + bool QuickHull::affineFilter(std::vector const& subset, uint_fast64_t const& item, std::vector> const& points){ + EigenMatrix vectorMatrix(vertices[item].dimension()+1, subset.size() + 1); + for (uint_fast64_t i = 0; i < subset.size(); ++i){ + vectorMatrix.col(i) << vertices[subset[i]], storm::utility::one(); + } + vectorMatrix.col(subset.size()) << vertices[item], storm::utility::one(); + return (vectorMatrix.fullPivLu().rank() > subset.size()); + } + + + template + bool QuickHull::findInitialVertices(std::vector& points, std::vector& verticesOfInitialPolytope, uint_fast64_t& minMaxVertices) const{ + const uint_fast64_t dimension = points[0].dimension(); + if(points.size() < dimension + 1){ + //not enough points to obtain a (non-degenerated) polytope + return false; + } + const uint_fast64_t candidatesToFind = std::min(2*dimension, points.size()); + uint_fast64_t candidatesFound = 0; + storm::storage::BitVector consideredPoints(points.size(), true); + while(candidatesFound < candidatesToFind && !consideredPoints.empty()) { + for(uint_fast64_t currDim=0; currDim points[pointIndex](currDim)){ + minIndex = pointIndex; + } + if(points[maxIndex](currDim) < points[pointIndex](currDim)){ + maxIndex = pointIndex; + } + } + consideredPoints.set(minIndex, false); + consideredPoints.set(maxIndex, false); + } + //Found candidates. Now swap them to the front. + consideredPoints = ~consideredPoints; + const uint_fast64_t newNumberOfCandidates = consideredPoints.getNumberOfSetBits(); + assert(newNumberOfCandidates > 0); + if(newNumberOfCandidates < points.size()){ + uint_fast64_t nextPointToMove = consideredPoints.getNextSetIndex(newNumberOfCandidates); + for(uint_fast64_t indexAtFront = candidatesFound; indexAtFront < newNumberOfCandidates; ++indexAtFront){ + if(!consideredPoints.get(indexAtFront)) { + assert(nextPointToMove != consideredPoints.size()); + std::swap(points[indexAtFront], points[nextPointToMove]); + nextPointToMove = consideredPoints.getNextSetIndex(nextPointToMove+1); + consideredPoints.set(indexAtFront); + } + } + assert(nextPointToMove == consideredPoints.size()); + } + if(candidatesFound==0){ + //We are in the first iteration. It holds that the first newNumberOfCandidates points will be vertices of the final polytope + minMaxVertices = newNumberOfCandidates; + } + candidatesFound = newNumberOfCandidates; + consideredPoints = ~consideredPoints; + } + + storm::storage::geometry::SubsetEnumerator>> subsetEnum(points.size(), dimension+1, points, affineFilter); + if(subsetEnum.setToFirstSubset()){ + verticesOfInitialPolytope = subsetEnum.getCurrentSubset(); + return true; + } else { + return false; + } + } + + template + std::vector::Facet> computeInitialFacets(std::vector const& points, std::vector const& verticesOfInitialPolytope, EigenVector const& insidePoint) const { + const uint_fast64_t dimension = points.front().rows(); + assert(verticesOfInitialPolytope.size() == dimension + 1); + std::vector result; + result.reserve(dimension + 1); + storm::storage::geometry::SubsetEnumerator<> subsetEnum(verticesOfInitialPolytope.size(), dimension); + if(!subsetEnum.setToFirstSubset()){ + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not find an initial subset."); + } + do{ + Facet newFacet; + // set the points that lie on the new facet + std::vector const& subset(subsetEnum.getCurrentSubset()); + newFacet.points.reserve(subset.size()); + for(uint_fast64_t i : subset){ + newFacet.points.push_back(verticesOfInitialPolytope[i]); + } + //neighbors: these are always the remaining facets + newFacet.neighbors.reserve(dimension); + for(uint_fast64_t i = 0; i < dimension+1; ++i){ + if(i != facets.size()){ //initFacets.size() will be the index of this new facet! + newFacet.neighbors.push_back(i); + } + } + // normal and offset: + computeNormalAndOffsetOfFacet(points, insidePoint, newFacet); + + facets.push_back(std::move(newFacet)); + } while(subsetEnum.incrementSubset()); + assert(facets.size()==dimension+1); + } + + template + void computeNormalAndOffsetOfFacet(std::vector const& points, EigenVector const& insidePoint, Facet& facet) const { + const uint_fast64_t dimension = points.front().rows() + assert(facet.points.size() == dimension); + EigenVector refPoint(facet.points.back()); + EigenMatrix constraints(dimension-1, dimension); + for(unsigned row = 0; row < dimension-1; ++row) { + constraints.row(row) = points[facet.points[row]] - refPoint; + } + facet.normal = constraints.fullPivLu().kernel(); + facet.offset = facet.normal.dot(refPoint); + + // invert the plane if the insidePoint is not contained in it + if(facet.normal.dot(insidePoint) > facet.offset) { + facet.normal = -facet.normal; + facet.offset = -facet.offset; + } + } + + + template + std::set QuickHull::getVisibleSet(std::vector const& facets, uint_fast64_t const& startIndex, EigenVector const& point) const { + std::set facetsToCheck; + std::set facetsChecked; + std::set visibleSet; + facetsChecked.insert(startIndex); + visibleSet.insert(startIndex); + for(uint_fast64_t i = 0; i < facets[startIndex].neighbors.size(); ++i){ + facetsToCheck.insert(facets[startIndex].neighbors[i]); + } + while (!facetsToCheck.empty()){ + auto elementIt = facetsToCheck.begin(); + if ( point.dot(facets[*elementIt].normal) > facets[*elementIt].offset) { + visibleSet.insert(*elementIt); + for(uint_fast64_t i = 0; i < facets[*elementIt].neighbors.size(); ++i){ + if (facetsChecked.find(facets[*elementIt].neighbors[i]) == facetsChecked.end()){ + facetsToCheck.insert(facets[*elementIt].neighbors[i]); + } + } + } + facetsChecked.insert(*elementIt); + facetsToCheck.erase(elementIt); + } + return visibleSet; + } + + template + void QuickHull::setNeighborhoodOfNewFacets(std::vector& facets, uint_fast64_t firstNewFacet, uint_fast64_t dimension) const{ + for(uint_fast64_t currentFacet = firstNewFacet; currentFacet < facets.size(); ++currentFacet){ + for(uint_fast64_t otherFacet = currentFacet + 1; otherFacet < facets.size(); ++otherFacet){ + std::vector commonPoints = getCommonPoints(facets[currentFacet], facets[otherFacet]); + if (commonPoints.size() >= dimension-1){ + facets[currentFacet].neighbors.push_back(otherFacet); + facets[otherFacet].neighbors.push_back(currentFacet); + } + } + } + } + + template + void QuickHull::replaceFacetNeighbor(std::vector& facets, uint_fast64_t oldFacetIndex, uint_fast64_t newFacetIndex, uint_fast64_t neighborIndex) const { + uint_fast64_t index = 0; + while(facets[neighborIndex].neighbors[index] != oldFacetIndex && index < facets[neighborIndex].neighbors.size()){ + ++index; + } + if (index < facets[neighborIndex].neighbors.size()){ + facets[neighborIndex].neighbors[index] = newFacetIndex; + } + } + + template + void QuickHull::computeOutsideSetOfFacet(Facet& facet, storm::storage::BitVector& currentOutsidePoints, std::vector const& points) const { + Number maxMultiplicationResult = facet.hyperplane.offset(); + for(uint_fast64_t pointIndex = currentOutsidePoints) { + ValueType multiplicationResult = points[pointIndex].dot(facet.normal); + if( multiplicationResult > facet.hyperplane.offset() ) { + currentOutsidePoints.set(pointIndex, false); // we already know that the point lies outside so it can be ignored for future facets + facet.outsideSet.push_back(pointIndex); + if (multiplicationResult > maxMultiplicationResult) { + maxMultiplicationResult = multiplicationResult; + facet.maxOutsidePointIndex = pointIndex; + } + } + } + } + + template + std::vector QuickHull::getCommonPoints(Facet const& lhs, Facet const& rhs) const{ + std::vector commonPoints; + for(uint_fast64_t refPoint = 0; refPoint < lhs.points.size(); ++refPoint){ + for(uint_fast64_t currentPoint = 0; currentPoint < rhs.points.size(); ++currentPoint){ + if (lhs.points[refPoint] == rhs.points[currentPoint]){ + commonPoints.push_back(lhs.points[refPoint]); + } + } + } + return commonPoints; + } + + template + std::set QuickHull::getInvisibleNeighbors( std::vector& facets, std::set const& visibleSet) const { + std::set invisibleNeighbors; + for(auto currentFacetIt = visibleSet.begin(); currentFacetIt != visibleSet.end(); ++currentFacetIt){ + for(uint_fast64_t currentNeighbor = 0; currentNeighbor < facets[*currentFacetIt].neighbors.size(); ++currentNeighbor){ + if (visibleSet.find(facets[*currentFacetIt].neighbors[currentNeighbor]) == visibleSet.end()){ + invisibleNeighbors.insert(facets[*currentFacetIt].neighbors[currentNeighbor]); + } + } + } + return invisibleNeighbors; + } + + template + void QuickHull::enlargeIncompleteResult(std::vector const& facets, storm::storage::BitVector const& currentFacets, std::vector const& points, bool generateRelevantVerticesAndVertexSets){ + PTERM_TRACE("Enlarging incomplete Result of QuickHull"); + //Obtain the set of outside points + std::vector outsidePoints; + for(uint_fast64_t facetIndex = currentFacets.find_first(); facetIndex != currentFacets.npos; facetIndex = currentFacets.find_next(facetIndex)){ + outsidePoints.insert(outsidePoints.end(), facets[facetIndex].outsideSet.begin(), facets[facetIndex].outsideSet.end()); + } + + //Now adjust the offsets of all the hyperplanes such that they contain each outside point + for(uint_fast64_t planeIndex=0; planeIndex < this->mResultMatrix.rows(); ++planeIndex){ + Number maxValue = this->mResultVector(planeIndex); + for (uint_fast64_t outsidePointIndex : outsidePoints){ + Number currValue = this->mResultMatrix.row(planeIndex) * (points[outsidePointIndex].rawCoordinates()); + maxValue = std::max(maxValue, currValue); + } + if (pterm::Settings::instance().getQuickHullRoundingMode() != pterm::Settings::QuickHullRoundingMode::NONE) { + maxValue = NumberReduction::instance().roundUp(maxValue); + } + this->mResultVector(planeIndex) = maxValue; + } + + if(generateRelevantVerticesAndVertexSets){ + /* Note: The adjustment of the offset might introduce redundant halfspaces. + * It is also not clear if it suffices to consider intersection points of hyperplanes that intersected in the original polytope + * + * Our solution is to convert the resulting h polytope back to a v poly + */ + + PTermHPolytope enlargedPoly(std::move(this->mResultMatrix), std::move(this->mResultVector)); + HyperplaneEnumeration he; + if(!he.generateVerticesFromHalfspaces(enlargedPoly, true)){ + PTERM_ERROR("Could not find the vertices of the enlarged Polytope"); + } + this->mResultMatrix = std::move(he.getRelevantMatrix()); + this->mResultVector = std::move(he.getRelevantVector()); + this->mRelevantVertices = std::move(he.getResultVertices()); + this->mVertexSets = std::move(he.getVertexSets()); + + } + } + } +} + +/* + * File: AbstractVtoHConverter.tpp + * Author: tim quatmann + * + * Created on Februrary 2, 2016, 1:06 PM + */ + +#include "AbstractVtoHConverter.h" + +namespace hypro { + namespace pterm{ + + template + EigenMatrix& AbstractVtoHConverter::getResultMatrix() { + return this->mResultMatrix; + } + + template + EigenVector& AbstractVtoHConverter::getResultVector() { + return this->mResultVector; + } + + template + std::vector& AbstractVtoHConverter::getRelevantVertices() { + return this->mRelevantVertices; + } + + template + std::vector>& AbstractVtoHConverter::getVertexSets() { + return this->mVertexSets; + } + + } +} + diff --git a/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h b/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h new file mode 100644 index 000000000..69359d1c9 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h @@ -0,0 +1,165 @@ +#ifndef STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_QUICKHULL_H_ +#define STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_QUICKHULL_H_ + +#include "storm/storage/BitVector.h" + +namespace storm { + namespace storage { + namespace geometry { + + template< typename ValueType> + class QuickHull { + public: + + typedef StormEigen::Matrix EigenMatrix; + typedef StormEigen::Matrix EigenVector; + + QuickHull() = default; + ~QuickHull() = default; + + + /* + * Generates the halfspaces of the given set of Points by the QuickHull-algorithm + * If the given flag is true, this method will also compute + * * the minimal set of vertices which represent the given polytope (can be used to remove redundant vertices), and + * * for each hyperplane, the set of (non-redundant) vertices that lie on each hyperplane. + * + * Use the provided getter methods to retrieve the results + * + * @return true iff conversion was successful. + */ + void generateHalfspacesFromPoints(std::vector const& points, bool generateRelevantVerticesAndVertexSets); + + + EigenMatrix& getResultMatrix(); + + EigenVector& getResultVector(); + + /*! + * Returns the set of vertices which are not redundant + * @note the returned vector is empty if the corresponding flag was false + */ + std::vector & getRelevantVertices(); + + /*! + * Returns for each hyperplane the set of vertices that lie on that hyperplane. + * A vertex is given as an index in the relevantVertices vector. + * @note the returned vector is empty if the corresponding flag was false + */ + std::vector>& getVertexSets(); + + + private: + + struct Facet { + EigenVector normal; + ValueType offset; + std::vector points; + std::vector neighbors; + // maxOutsidePointIndex and outsideSet will be set in Quickhull algorithm + std::vector outsideSet; + uint_fast64_t maxOutsidePointIndex; + }; + + /* + * Returns true if the vertices with index of subset and item are affine independent + * Note that this filter also works for dimension()+1 many vertices + */ + static bool affineFilter(std::vector const& subset, uint_fast64_t const& item, std::vector const& vertices); + + + /*! + * finds a set of vertices that correspond to a (hopefully) large V polytope. + * + * @param points The set of points from which vertices are picked. Note that the order of the points might be changed when calling this!! + * @param verticesOfInitialPolytope Will be set to the resulting vertices (represented by indices w.r.t. the given points) + * @param minMaxVertices after calling this, the first 'minMaxVertices' points will have an extreme value in at least one dimension + * @return true if the method was successful. False if the given points are affine dependend, i.e. the polytope is degenerated. + */ + bool findInitialVertices(std::vector& points, std::vector& verticesOfInitialPolytope, uint_fast64_t& minMaxVertices) const; + + /*! + * Computes the initial facets out of the given dimension+1 initial vertices + */ + std::vector computeInitialFacets(std::vector const& points, std::vector const& verticesOfInitialPolytope, EigenVector const& insidePoint) const; + + + // Computes the normal vector and the offset of the given facet from the (dimension many) points specified in the facet. + // The insidePoint specifies the orientation of the facet. + void computeNormalAndOffsetOfFacet(std::vector const& points, EigenVector const& insidePoint, Facet& facet) const; + + /* + * Extends the given mesh using the QuickHull-algorithm + * For optimization reasons a point thats inside of the initial polytope but on none of the facets has to be provided. + + * @return true iff all consideredPoints are now contained by the mesh + */ + void extendMesh(std::vector& points, + storm::storage::BitVector& consideredPoints, + std::vector>& facets, + storm::storage::BitVector& currentFacets, + vector_t& insidePoint, + uint_fast64_t& currentNumOfVertices) const; + + /*! + * Uses the provided mesh to generate a HPolytope (in form of a matrix and a vector) + * If the given flag is true, this method will also compute + * * the minimal set of vertices which represent the given vPoly (can be used to remove redundant vertices), and + * * for each hyperplane, the set of (non-redundant) vertices that lie on each hyperplane. + * + */ + void getPolytopeFromMesh(std::vector const& points, std::vector> const& facets, storm::storage::BitVector const& currentFacets, bool generateRelevantVerticesAndVertexSets); + + + /* + * Returns the set of facets visible from point starting with the facet with index startIndex and recursively testing all neighbors + */ + std::set getVisibleSet(std::vector> const& facets, uint_fast64_t const& startIndex, EigenVector const& point) const; + + /* + * Sets neighborhood for all facets with index >= firstNewFacet in facets + */ + void setNeighborhoodOfNewFacets(std::vector>& facets, uint_fast64_t firstNewFacet, uint_fast64_t dimension) const; + + /* + * replaces oldFacet by newFacet in the neighborhood of neighbor + */ + void replaceFacetNeighbor(std::vector>& facets, uint_fast64_t oldFacetIndex, uint_fast64_t newFacetIndex, uint_fast64_t neighborIndex) const; + + /* + * computes the outside set of the given facet + */ + void computeOutsideSetOfFacet(Facet& facet, storm::storage::BitVector& currentOutsidePoints, std::vector const& points) const; + + /* + * returns common points of lhs and rhs + */ + std::vector getCommonPoints(Facet const& lhs, Facet const& rhs) const; + + /* + * computes all neighbors that are not in the visibleSet + */ + std::set getInvisibleNeighbors( std::vector>& facets, std::set const& visibleSet) const; + + bool facetHasNeighborWithIndex(Facet const& facet, uint_fast64_t searchIndex) const; + + /* + * Enlarges the result such that all points that are still outside will be contained in the resulting polytope. + */ + void enlargeIncompleteResult(std::vector> const& facets, storm::storage::BitVector const& currentFacets, std::vector const& points, bool generateRelevantVerticesAndVertexSets); + + + hypro::matrix_t resultMatrix; + hypro::vector_t resultVector; + std::vector relevantVertices; + std::vector > vertexSets; + + }; + } +} + + + + + +#endif STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_QUICKHULL_H_ \ No newline at end of file diff --git a/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp b/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp new file mode 100644 index 000000000..1ba9e91c7 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp @@ -0,0 +1,93 @@ +#include "SubsetEnumerator.h" + + +#include "storm/utility/eigen.h" + +namespace storm { + namespace storage { + namespace geometry { + + template< typename DataType> + SubsetEnumerator::SubsetEnumerator(uint_fast64_t n, uint_fast64_t k, DataType const& data, SubsetFilter subsetFilter) : n(n), k(k), data(data), filter(subsetFilter){ + //Intentionally left empty + } + + template< typename DataType> + SubsetEnumerator::~SubsetEnumerator(){ + //Intentionally left empty + } + + template< typename DataType> + std::vector const& SubsetEnumerator::getCurrentSubset(){ + return this->current; + } + + template< typename DataType> + bool SubsetEnumerator::setToFirstSubset(){ + if(n < k) return false; + //set the upper boundaries first. + upperBoundaries.clear(); + upperBoundaries.reserve(k); + for(uint_fast64_t bound = (n - k); bound < n; ++bound){ + upperBoundaries.push_back(bound); + } + //now set the current subset to the very first one + current.clear(); + current.reserve(k); + uint_fast64_t newItem = 0; + while(current.size()!= k && newItem<=upperBoundaries[current.size()]){ + //Check if it is okay to add the new item... + if(filter(current, newItem, data)){ + current.push_back(newItem); + } + ++newItem; + } + //Note that we only insert things into the vector if it is "okay" to do so. + //Hence, we have failed iff we were not able to insert k elements. + return current.size() == k; + } + + template< typename DataType> + bool SubsetEnumerator::incrementSubset(){ + //The currentSelection will be the numbers that are already inside of our new subset. + std::vector currentSelection(current); + currentSelection.pop_back(); + uint_fast64_t pos = k-1; + while (true){ + //check whether we can increment at the current position + if(current[pos] == upperBoundaries[pos]){ + if(pos==0){ + //we already moved to the very left and still can not increment. + //Hence, we are already at the last subset + return false; + } + currentSelection.pop_back(); + --pos; + } else { + ++current[pos]; + //check if the new subset is inside our filter + if(filter(currentSelection, current[pos], data)){ + //it is, so add it and go on with the position on the right + currentSelection.push_back(current[pos]); + ++pos; + if(pos == k){ + //we are already at the very right. + //Hence, we have found our new subset of size k + return true; + } + //initialize the value at the new position + current[pos] = current[pos-1]; + } + } + } + } + + template< typename DataType> + bool SubsetEnumerator::trueFilter(std::vector const& subset, uint_fast64_t const& item, DataType const& data){ + return true; + } + + template class SubsetEnumerator; + } + } +} \ No newline at end of file diff --git a/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h b/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h new file mode 100644 index 000000000..2d30d9ce3 --- /dev/null +++ b/src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h @@ -0,0 +1,69 @@ +#ifndef STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_SUBSETENUMERATOR_H_ +#define STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_SUBSETENUMERATOR_H_ + +#include + + +namespace storm { + namespace storage { + namespace geometry { + /*! + * This class can be used to enumerate all k-sized subsets of {0,...,n-1}. + * A subset is represented as a vector of ascending numbers. + * Example: (n=5, k=3) + * [0,1,2] --> [0,1,3] --> [0,1,4] --> [0,2,3] --> [0,2,4] --> [0,3,4] --> + * [1,2,3] --> [1,2,4] --> [1,3,4] --> [2,3,4] + * A filter can be given which should answer true iff it is ok to add a given + * item to a given subset. + * Example: (n=5, k=3, filter answers false iff 0 and 2 would be in the + * resulting subset): + * [0,1,3] --> [0,1,4] --> [0,3,4] --> [1,2,3] --> [1,2,4] --> [1,3,4] --> [2,3,4] + */ + template + class SubsetEnumerator { + + public: + + //A typedef for the filter. + //Note that the function will be called with subset.size() in {0, ..., k-1}. + typedef bool (*SubsetFilter)(std::vector const& subset, uint_fast64_t const& item, DataType const& data); + + /* + * Constructs a subset enumerator that can enumerate all k-sized Subsets of {0,...,n-1} + * The given filter can be used to skip certain subsets. + * @note call "setToFirstSubset()" before retrieving the first subset + */ + SubsetEnumerator(uint_fast64_t n, uint_fast64_t k, DataType const& data = DataType(), SubsetFilter subsetFilter = trueFilter); + + ~SubsetEnumerator(); + + //returns the current subset of size k. + //Arbitrary behavior if setToFirstSubset or incrementSubset returned false or have never been executed + std::vector const& getCurrentSubset(); + + // Sets the current subset to the very first one. + // @note Needs to be called initially. + // Returns true iff there actually is a first subset and false if not (e.g. when n const& subset, uint_fast64_t const& item, DataType const& data); + + private: + uint_fast64_t n; //the size of the source set + uint_fast64_t k; //the size of the desired subsets + DataType const &data; //The data which is given as additional information when invoking the filter + SubsetFilter filter; //returns true iff it is okay to insert a new element + std::vector current; //the current subset + std::vector upperBoundaries; //will always be [n-k, ..., n-1]. Used to easily check whether we can increment the subset at a given position + }; + } + } +} + +#endif /* STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_SUBSETENUMERATOR_H_ */ \ No newline at end of file