started on native polytopes

src/storm/adapters/EigenAdapter.cpp

@@ -19,12 +19,29 @@ namespace storm {
             result->setFromTriplets(triplets.begin(), triplets.end());
             return result;
         }
-        
+
+        template <typename ValueType>
+        std::vector<ValueType> EigenAdapter::toStdVector(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> const& v){
+            return std::vector<ValueType>(v.data(), v.data() + v.rows());
+        }
+
+        template <typename ValueType>
+        StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<ValueType> const& v){
+            return StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(v.data(), v.size());
+        }
+
         template std::unique_ptr<StormEigen::SparseMatrix<double>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<double> const& matrix);
+        template std::vector<double> EigenAdapter::toStdVector(StormEigen::Matrix<double, StormEigen::Dynamic, 1> const& v);
+        template StormEigen::Matrix<double, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<double> const& v);
 
 #ifdef STORM_HAVE_CARL
         template std::unique_ptr<StormEigen::SparseMatrix<storm::RationalNumber>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix);
+        template std::vector<storm::RationalNumber> EigenAdapter::toStdVector(StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1> const& v);
+        template StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalNumber> const& v);
+
         template std::unique_ptr<StormEigen::SparseMatrix<storm::RationalFunction>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix);
+        template std::vector<storm::RationalFunction> EigenAdapter::toStdVector(StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1> const& v);
+        template StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalFunction> const& v);
 #endif
     }
 }

src/storm/adapters/EigenAdapter.h

@@ -17,6 +17,12 @@ namespace storm {
              */
             template<class ValueType>
             static std::unique_ptr<StormEigen::SparseMatrix<ValueType>> toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix);
+
+            template <typename ValueType>
+            static std::vector<ValueType> toStdVector(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> const& v);
+
+            template <typename ValueType>
+            static StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> toEigenVector(std::vector<ValueType> const& v);
         };
         
     }

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 <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(std::vector<Halfspace<ValueType>> 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<ValueType>::toEigenVector(halfspaces[row].normal());
+                    }
+                    emptyStatus = EmptyStatus::Unknown;
+                }
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(std::vector<Point> const& points) {
+                if(points.empty()){
+                    emptyStatus = EmptyStatus::Empty;
+                } else {
+                    std::vector<EigenVector> eigenPoints;
+                    eigenPoints.reserve(points.size());
+                    for(auto const& p : points){
+                        eigenPoints.emplace_back(storm::adapters::EigenAdapter<ValueType>::toEigenVector(p));
+                    }
+
+                    // todo: quickhull(eigenPoints)
+
+                    emptyStatus = EmptyStatus::Nonempty;
+                }
+            }
+            
+            template <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<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 NativePolytope where halfspaces AND points are given. The points will be ignored.");
+                    return std::make_shared<NativePolytope<ValueType>>(*halfspaces);
+                } else if(points) {
+                    return std::make_shared<NativePolytope<ValueType>>(*points);
+                }
+                STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Creating a NativePolytope but no representation was given.");
+                return nullptr;
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(NativePolytope<ValueType> const& other) : emptyStatus(other.emptyStatus), A(other.A), b(other.b) {
+                // Intentionally left empty
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(NativePolytope<ValueType>&& other) : emptyStatus(std::move(other.emptyStatus)), A(std::move(other.A)), b(std::move(other.b)) {
+                // Intentionally left empty
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(EmptyStatus const& emptyStatus, EigenMatrix const& halfspaceMatrix, EigenVector const& halfspaceVector) : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) {
+                // Intentionally left empty
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::NativePolytope(EmptyStatus&& emptyStatus, EigenMatrix&& halfspaceMatrix, EigenVector&& halfspaceVector) : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) {
+                // Intentionally left empty
+            }
+            
+            template <typename ValueType>
+            NativePolytope<ValueType>::~NativePolytope() {
+                // Intentionally left empty
+            }
+            
+            template <typename ValueType>
+            std::vector<typename Polytope<ValueType>::Point> NativePolytope<ValueType>::getVertices() const {
+                std::vector<hypro::Point<ValueType>> eigenVertices = getEigenVertices();
+                std::vector<Point> result;
+                result.reserve(eigenVertices.size());
+                for(auto const& p : eigenVertices) {
+                    result.push_back(storm::adapters::EigenAdapter<ValueType>::toStdVector(p));
+                }
+                return result;
+            }
+            
+            template <typename ValueType>
+            std::vector<Halfspace<ValueType>> NativePolytope<ValueType>::getHalfspaces() const {
+                std::vector<Halfspace<ValueType>> result;
+                result.reserve(A.rows());
+
+                for(uint_fast64_t row=0; row < A.rows(); ++row){
+                    result.emplace_back(storm::adapters::EigenAdapter<ValueType>::toStdVector(A.row(row)), b.row(row));
+                }
+            }
+
+            template <typename ValueType>
+            bool NativePolytope<ValueType>::isEmpty() const {
+                if(emptyStatus == emptyStatus::Unknown) {
+                    storm::expressions::ExpressionManager manager;
+                    std::unique_ptr<storm::solver::SmtSolver> 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 <typename ValueType>
+            bool NativePolytope<ValueType>::isUniversal() const {
+                return A.rows()==0;
+            }
+            
+            template <typename ValueType>
+            bool NativePolytope<ValueType>::contains(Point const& point) const{
+                EigenVector x = storm::adapters::EigenAdapter<ValueType>::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 <typename ValueType>
+            bool NativePolytope<ValueType>::contains(std::shared_ptr<Polytope<ValueType>> 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<storm::solver::SmtSolver> 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<NativePolytope<ValueType> 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 <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::intersection(std::shared_ptr<Polytope<ValueType>> 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<ValueType> const& nativeRhs = dynamic_cast<NativePolytope<ValueType> 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<NativePolytope<ValueType>>(std::move(resultA), std::move(resultb));
+            }
+            
+            template <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::intersection(Halfspace<ValueType> const& halfspace) const{
+                EigenMatrix resultA(A.rows() + 1, A.cols());
+                resultA << A,
+                           storm::adapters::EigenAdapter<ValueType>::toEigenVector(halfspace.normal());
+                EigenVector resultb(resultA.rows());
+                resultb << b,
+                           halfspace.offset();
+                return std::make_shared<NativePolytope<ValueType>>(std::move(resultA), std::move(resultb));
+            }
+            
+            template <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::convexUnion(std::shared_ptr<Polytope<ValueType>> 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<NativePolytope<ValueType>>(dynamic_cast<NativePolytope<ValueType> const&>(*rhs));
+                } else if (rhs->isEmpty()) {
+                    return std::make_shared<NativePolytope<ValueType>>(*this);
+                }
+                if(this->isUniversal() || rhs->isUniversal) {
+                    return std::make_shared<NativePolytope<ValueType>>(*this);
+                }
+
+                STORM_LOG_WARN("Implementation of convex union of two polytopes only works if the polytopes are bounded.");
+
+                std::vector<EigenVector> rhsVertices = dynamic_cast<NativePolytope<ValueType> const&>(*rhs).getEigenVertices();
+                std::vector<EigenVector> 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 <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::minkowskiSum(std::shared_ptr<Polytope<ValueType>> 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<ValueType> const& nativeRhs = dynamic_cast<NativePolytope<ValueType> const&>(*rhs);
+
+                return std::make_shared<NativePolytope<ValueType>>(internPolytope.minkowskiSum(dynamic_cast<NativePolytope<ValueType> const&>(*rhs).internPolytope));
+            }
+            
+            template <typename ValueType>
+            std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::affineTransformation(std::vector<Point> 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<ValueType> 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<NativePolytope<ValueType>>(internPolytope.affineTransformation(std::move(hyproMatrix), storm::adapters::toNative(vector)));
+            }
+            
+            template <typename ValueType>
+            std::pair<typename NativePolytope<ValueType>::Point, bool> NativePolytope<ValueType>::optimize(Point const& direction) const {
+                hypro::EvaluationResult<ValueType> 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 <typename ValueType>
+            bool NativePolytope<ValueType>::isNativePolytope() const {
+                return true;
+            }
+            template <typename ValueType>
+            std::vector<typename NativePolytope::ValueType>::EigenVector> NativePolytope<ValueType>::getEigenVertices() const {
+                // todo: invoke conversion
+                return std::vector<EigenVector>();
+            }
+
+            template <typename ValueType>
+            std::vector<storm::expressions::Variable> NativePolytope<ValueType>::declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const {
+                std::vector<storm::expressions::Variable> 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 <typename ValueType>
+            storm::expressions::Expression NativePolytope<ValueType>::getConstraints(storm::expressions::ExpressionManager& manager, std::vector<storm::expressions::Variable> 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<double>;
+#ifdef STORM_HAVE_CARL
+            template class NativePolytope<storm::RationalNumber>;
+#endif
+        }
+    }
+}
+#endif

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 <typename ValueType>
+            class NativePolytope : public Polytope<ValueType> {
+            public:
+                
+                typedef typename Polytope<ValueType>::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<Polytope<ValueType>> create(boost::optional<std::vector<Halfspace<ValueType>>> const& halfspaces,
+                                                                   boost::optional<std::vector<Point>> const& points);
+                
+                /*!
+                 * Creates a NativePolytope from the given halfspaces
+                 * The resulting polytope is defined as the intersection of the halfspaces.
+                 */
+                NativePolytope(std::vector<Halfspace<ValueType>> const& halfspaces);
+                
+                /*!
+                 * Creates a NativePolytope from the given points.
+                 * The resulting polytope is defined as the convex hull of the points'
+                 */
+                NativePolytope(std::vector<Point> const& points);
+                
+                /*!
+                 * Copy and move constructors
+                 */
+                NativePolytope(NativePolytope<ValueType> const& other);
+                NativePolytope(NativePolytope<ValueType>&& 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<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 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<Polytope<ValueType>> affineTransformation(std::vector<Point> 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<Point, bool> optimize(Point const& direction) const override;
+            
+                virtual bool isNativePolytope() const override;
+                
+            private:
+                typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
+                typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenVector;
+
+                // returns the vertices of this polytope as EigenVectors
+                std::vector<EigenVector> getEigenVertices() const;
+
+                // declares one variable for each constraint and returns the obtained variables.
+                std::vector<storm::expressions::Variable> 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<storm::expressions::Variable> 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_ */

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<Number>::insert(hypro::Hyperplane<Number>const & hyperplane, std::vector<std::size_t>const* indexList) {
+            return this->insert(hyperplane.normal(), hyperplane.offset(), indexList);
+        }
+        
+        template< typename Number>
+        bool HyperplaneCollector<Number>::insert(hypro::vector_t<Number> const& normal, Number const& offset, std::vector<std::size_t>const* indexList) {
+            hypro::vector_t<Number> copyNormal(normal);
+            Number copyOffset(offset);
+            return this->insert(std::move(copyNormal), std::move(copyOffset), indexList);
+        }
+        
+        template< typename Number>
+        bool HyperplaneCollector<Number>::insert(hypro::vector_t<Number> && normal, Number && offset, std::vector<std::size_t>const* indexList) {
+            //Normalize
+            Number infinityNorm = normal.template lpNorm<Eigen::Infinity>();
+            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<std::size_t>())).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::matrix_t<Number>, hypro::vector_t<Number>> HyperplaneCollector<Number>::getCollectedHyperplanesAsMatrixVector() const{
+            assert(!this->mMap.empty());
+            hypro::matrix_t<Number> A(this->mMap.size(), this->mMap.begin()->first.first.rows());
+            hypro::vector_t<Number> 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::matrix_t<Number>, hypro::vector_t<Number>>(std::move(A), std::move(b)); 
+        }
+
+        template< typename Number>
+        std::vector<std::vector<std::size_t>> HyperplaneCollector<Number>::getIndexLists() const{
+            assert(!this->mMap.empty());
+            std::vector<std::vector<std::size_t>> 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<Number>::numOfCollectedHyperplanes() const {
+            return this->mMap.size();
+        }
+    }
+}
+

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 <unordered_map>
+#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<Number> const& hyperplane, std::vector<std::size_t>const* indexList = nullptr);
+            bool insert(hypro::vector_t<Number> const& normal, Number const& offset, std::vector<std::size_t>const* indexList = nullptr);
+            bool insert(hypro::vector_t<Number>&& normal, Number && offset, std::vector<std::size_t>const* indexList = nullptr);
+            
+            std::pair<hypro::matrix_t<Number>, hypro::vector_t<Number>> getCollectedHyperplanesAsMatrixVector() const;
+            //Note that the returned lists might contain dublicates.
+            std::vector<std::vector<std::size_t>> getIndexLists() const;
+            
+            std::size_t numOfCollectedHyperplanes() const;
+            
+        private:
+            typedef std::pair<hypro::vector_t<Number>, Number> NormalOffset; 
+            class NormalOffsetHash{
+                public:
+                    std::size_t operator()(NormalOffset const& ns) const {
+                        std::size_t seed = std::hash<hypro::vector_t<Number>>()(ns.first);
+                        carl::hash_add(seed, std::hash<Number>()(ns.second));
+                        return seed;
+                    }
+            };
+            typedef typename std::unordered_map<NormalOffset, std::vector<std::size_t>, NormalOffsetHash>::key_type KeyType;
+            typedef typename std::unordered_map<NormalOffset, std::vector<std::size_t>, NormalOffsetHash>::value_type ValueType;
+            
+            std::unordered_map<NormalOffset, std::vector<std::size_t>, NormalOffsetHash> mMap;
+            
+            
+        };
+    }
+}
+
+
+#include "HyperplaneCollector.tpp"

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<typename Number>
+        bool HyperplaneEnumeration<Number>::generateVerticesFromHalfspaces(PTermHPolytope<Number> const& hPoly, bool generateRelevantHyperplanesAndVertexSets){
+            PTERM_DEBUG("Invoked generateVerticesFromHalfspaces with " << hPoly.getMatrix().rows() << " hyperplanes. Dimension is " << hPoly.dimension());
+            std::unordered_map<Point<Number>, std::set<std::size_t>> vertexCollector;
+            hypro::pterm::SubsetEnumerator<hypro::matrix_t<Number>> subsetEnum(hPoly.getMatrix().rows(), hPoly.dimension(), hPoly.getMatrix(), this->linearDependenciesFilter);
+            if(subsetEnum.setToFirstSubset()){
+                do{
+                    std::vector<std::size_t> const& subset = subsetEnum.getCurrentSubset();
+                    std::pair<hypro::matrix_t<Number>, hypro::vector_t<Number>> subHPolytope(hPoly.getSubHPolytope(subset));
+                    Point<Number> point(subHPolytope.first.fullPivLu().solve(subHPolytope.second));
+                    //Point<Number> 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<Point<Number>, std::set<std::size_t>>::value_type(std::move(point), std::set<std::size_t>())).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<std::size_t> 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<Number> hyperplaneCollector;
+                for(std::size_t hyperplaneIndex = 0; hyperplaneIndex < verticesOnHyperplaneCounter.size(); ++hyperplaneIndex){
+                    if(verticesOnHyperplaneCounter[hyperplaneIndex] >= hPoly.dimension()){
+                        std::vector<std::size_t> 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<std::size_t> oldToNewIndexMapping (hPoly.getMatrix().rows(), hPoly.getMatrix().rows()); //Initialize with some illegal value
+                std::vector<std::vector<std::size_t>> 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<std::set<std::size_t>> 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<Number>();
+                this->mRelevantVector = hypro::vector_t<Number>();
+            }
+            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 <typename Number>
+    bool HyperplaneEnumeration<Number>::linearDependenciesFilter(std::vector<std::size_t> const& subset, std::size_t const& item, hypro::matrix_t<Number> const& A) {
+        hypro::matrix_t<Number> 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<matrix_t<Number>> 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<typename Number>
+        std::vector<Point<Number>>& HyperplaneEnumeration<Number>::getResultVertices() {
+            return this->mResultVertices;
+        }
+        
+        template<typename Number>
+        hypro::matrix_t<Number>& HyperplaneEnumeration<Number>::getRelevantMatrix() {
+            return this->mRelevantMatrix;
+        }
+        
+        template<typename Number>
+        hypro::vector_t<Number>& HyperplaneEnumeration<Number>::getRelevantVector() {
+            return this->mRelevantVector;
+        }
+        
+        template<typename Number>
+        std::vector<std::vector<std::size_t>>& HyperplaneEnumeration<Number>::getVertexSets() {
+            return this->mVertexSets;
+        }
+        
+    }
+}
+

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<Number> const& hPoly, bool generateRelevantHyperplanesAndVertexSets);
+            
+            std::vector<Point<Number>>& 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<Number>& getRelevantMatrix();
+            hypro::vector_t<Number>& 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<std::vector<std::size_t>>& 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<std::size_t> const& subset, std::size_t const& item, hypro::matrix_t<Number> const& A);
+
+        private:
+            std::vector<Point<Number>> mResultVertices;
+            hypro::matrix_t<Number> mRelevantMatrix;
+            hypro::vector_t<Number> mRelevantVector;
+            std::vector<std::vector<std::size_t>> mVertexSets;
+        };
+    }
+}
+
+
+#include "HyperplaneEnumeration.tpp"

src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp

@@ -0,0 +1,505 @@
+#include "storm/storage/geometry/nativepolytopeconversion/QuickHull.h"
+
+#include <algorithm>
+
+#include "storm/utility/macros.h"
+
+#include "storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.h"
+#include "storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h"
+
+namespace hypro {
+    namespace pterm{
+        
+        template<typename ValueType>
+        void QuickHull<Number>::generateHalfspacesFromVertices(std::vector<EigenVector> 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<uint_fast64_t> 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<ValueType>(vertexIndices.size());
+
+            // Create the initial facets from the found vertices.
+            std::vector<Facet> 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<typename ValueType>
+        void QuickHull<Number>::extendMesh(std::vector<EigenVector>& points,
+                storm::storage::BitVector& consideredPoints,
+                std::vector<Facet>& facets,
+                storm::storage::BitVector& currentFacets,
+                vector_t<Number>& 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<uint_fast64_t> visibleSet = getVisibleSet(facets, facetCount, points[farAwayPointIndex]);
+                    std::set<uint_fast64_t> 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<EigenVector> 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<Number>(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<typename ValueType>
+        void QuickHull<Number>::getPolytopeFromMesh(std::vector<EigenVector> const& points, std::vector<Facet> const& facets, storm::storage::BitVector const& currentFacets, bool generateRelevantVerticesAndVertexSets){
+            
+            hypro::pterm::HyperplaneCollector<Number> 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<uint_fast64_t> hyperplanesOnVertexCounter(points.size(), 0);
+                for(auto& vertexVector : this->mVertexSets){
+                    std::set<uint_fast64_t> 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<hypro::EigenVector, std::vector<uint_fast64_t>> 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<hypro::EigenVector, std::vector<uint_fast64_t>>::value_type(points[vertexIndex], std::vector<uint_fast64_t>())).first;
+                        mapEntry->second.push_back(vertexIndex);
+                    }
+                }
+                //Fill in the relevant vertices and create a translation map from old to new indices
+                std::vector<uint_fast64_t> 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<uint_fast64_t> 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 <typename ValueType>
+        bool QuickHull<ValueType>::affineFilter(std::vector<uint_fast64_t> const& subset, uint_fast64_t const& item, std::vector<EigenVector<ValueType>> 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<ValueType>();
+            }
+            vectorMatrix.col(subset.size()) << vertices[item], storm::utility::one<ValueType>();
+            return (vectorMatrix.fullPivLu().rank() > subset.size());
+        }
+            
+        
+        template<typename ValueType>
+        bool QuickHull<Number>::findInitialVertices(std::vector<hypro::EigenVector>& points, std::vector<uint_fast64_t>& 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<dimension; ++currDim) {
+                    uint_fast64_t minIndex = *consideredPoints.begin();
+                    uint_fast64_t maxIndex = minIndex;
+                    for(uint_fast64_t pointIndex : consideredPoints){
+                        //Check if the current point is a new minimum or maximum at the current dimension
+                        if(points[minIndex](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<std::vector<EigenVector<ValueType>>> subsetEnum(points.size(), dimension+1, points, affineFilter);
+            if(subsetEnum.setToFirstSubset()){
+                verticesOfInitialPolytope = subsetEnum.getCurrentSubset();
+                return true;
+            } else {
+                return false;
+            }
+        }
+
+        template<typename ValueType>
+        std::vector<typename QuickHull<ValueType>::Facet> computeInitialFacets(std::vector<EigenVector> const& points, std::vector<uint_fast64_t> const& verticesOfInitialPolytope, EigenVector const& insidePoint) const {
+            const uint_fast64_t dimension = points.front().rows();
+            assert(verticesOfInitialPolytope.size() == dimension + 1);
+            std::vector<Facet> 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<uint_fast64_t> 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<typename ValueType>
+        void computeNormalAndOffsetOfFacet(std::vector<EigenVector> 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<typename ValueType>
+        std::set<uint_fast64_t> QuickHull<Number>::getVisibleSet(std::vector<Facet> const& facets, uint_fast64_t const& startIndex, EigenVector const& point) const {
+            std::set<uint_fast64_t> facetsToCheck;
+            std::set<uint_fast64_t> facetsChecked;
+            std::set<uint_fast64_t> 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<typename ValueType>
+        void QuickHull<Number>::setNeighborhoodOfNewFacets(std::vector<Facet>& 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<uint_fast64_t> commonPoints = getCommonPoints(facets[currentFacet], facets[otherFacet]); 
+                    if (commonPoints.size() >= dimension-1){
+                        facets[currentFacet].neighbors.push_back(otherFacet);
+                        facets[otherFacet].neighbors.push_back(currentFacet);
+                    }
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        void QuickHull<Number>::replaceFacetNeighbor(std::vector<Facet>& 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<typename ValueType>
+        void QuickHull<Number>::computeOutsideSetOfFacet(Facet& facet, storm::storage::BitVector& currentOutsidePoints, std::vector<EigenVector> 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<typename ValueType>
+        std::vector<uint_fast64_t> QuickHull<Number>::getCommonPoints(Facet const& lhs, Facet const& rhs) const{
+            std::vector<uint_fast64_t> 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<typename ValueType>
+        std::set<uint_fast64_t> QuickHull<Number>::getInvisibleNeighbors( std::vector<Facet>& facets, std::set<uint_fast64_t> const& visibleSet) const {
+            std::set<uint_fast64_t> 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<typename ValueType>
+        void QuickHull<Number>::enlargeIncompleteResult(std::vector<Facet> const& facets, storm::storage::BitVector const& currentFacets, std::vector<EigenVector> const& points,  bool generateRelevantVerticesAndVertexSets){
+            PTERM_TRACE("Enlarging incomplete Result of QuickHull");
+            //Obtain the set of outside points
+            std::vector<uint_fast64_t> 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<Number>::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<Number> enlargedPoly(std::move(this->mResultMatrix), std::move(this->mResultVector));
+                HyperplaneEnumeration<Number> 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<typename ValueType>
+        EigenMatrix& AbstractVtoHConverter<Number>::getResultMatrix() {
+            return this->mResultMatrix;
+        }
+        
+        template<typename ValueType>
+        EigenVector& AbstractVtoHConverter<Number>::getResultVector() {
+            return this->mResultVector;
+        }
+        
+        template<typename ValueType>
+        std::vector<EigenVector>& AbstractVtoHConverter<Number>::getRelevantVertices() {
+            return this->mRelevantVertices;
+        }
+        
+        template<typename ValueType>
+        std::vector<std::vector<uint_fast64_t>>& AbstractVtoHConverter<Number>::getVertexSets() {
+            return this->mVertexSets;
+        }
+        
+    }
+}
+

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<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
+            typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> 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<EigenVector> 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 <EigenVector>& 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<std::vector<std::uint_fast64_t>>& getVertexSets();
+
+                
+        private:
+
+            struct Facet {
+                EigenVector normal;
+                ValueType offset;
+                std::vector<uint_fast64_t> points;
+                std::vector<uint_fast64_t> neighbors;
+                // maxOutsidePointIndex and outsideSet will be set in Quickhull algorithm
+                std::vector<uint_fast64_t> 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<uint_fast64_t> const& subset, uint_fast64_t const& item, std::vector<EigenVector> 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<EigenVector>& points, std::vector<uint_fast64_t>& verticesOfInitialPolytope, uint_fast64_t& minMaxVertices) const;
+            
+            /*!
+             * Computes the initial facets out of the given dimension+1 initial vertices
+             */
+            std::vector<Facet> computeInitialFacets(std::vector<EigenVector> const& points, std::vector<uint_fast64_t> 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<EigenVector> 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<EigenVector>& points,
+                    storm::storage::BitVector& consideredPoints,
+                    std::vector<Facet<Number>>& facets,
+                    storm::storage::BitVector& currentFacets,
+                    vector_t<Number>& 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<EigenVector> const& points, std::vector<Facet<Number>> 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<uint_fast64_t> getVisibleSet(std::vector<Facet<Number>> 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<Facet<Number>>& facets, uint_fast64_t firstNewFacet, uint_fast64_t dimension) const;
+ 
+            /*
+             * replaces oldFacet by newFacet in the neighborhood of neighbor
+             */
+            void replaceFacetNeighbor(std::vector<Facet<Number>>& 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<Number>& facet, storm::storage::BitVector& currentOutsidePoints, std::vector<EigenVector> const& points) const;
+
+            /*
+             * returns common points of lhs and rhs
+             */
+            std::vector<uint_fast64_t> getCommonPoints(Facet<Number> const& lhs, Facet<Number> const& rhs) const;      
+ 
+            /*
+             * computes all neighbors that are not in the visibleSet
+             */
+            std::set<uint_fast64_t> getInvisibleNeighbors( std::vector<Facet<Number>>& facets, std::set<uint_fast64_t> const& visibleSet) const;
+
+            bool facetHasNeighborWithIndex(Facet<Number> 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<Facet<Number>> const& facets, storm::storage::BitVector const& currentFacets, std::vector<EigenVector> const& points,  bool generateRelevantVerticesAndVertexSets);
+
+
+            hypro::matrix_t<Number> resultMatrix;
+            hypro::vector_t<Number> resultVector;
+            std::vector <EigenVector> relevantVertices;
+            std::vector <std::vector<uint_fast64_t>> vertexSets;
+        
+        };
+    }
+}
+
+
+
+
+
+#endif STORM_STORAGE_GEOMETRY_NATIVEPOLYTOPECONVERSION_QUICKHULL_H_

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<DataType>::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<DataType>::~SubsetEnumerator(){
+                //Intentionally left empty
+            }
+
+            template< typename DataType>
+            std::vector<uint_fast64_t> const& SubsetEnumerator<DataType>::getCurrentSubset(){
+                return this->current;
+            }
+
+            template< typename DataType>
+            bool SubsetEnumerator<DataType>::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<DataType>::incrementSubset(){
+                //The currentSelection will be the numbers that are already inside of our new subset.
+                std::vector<uint_fast64_t> 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<DataType>::trueFilter(std::vector<uint_fast64_t> const& subset, uint_fast64_t const& item, DataType const& data){
+                return true;
+            }
+
+            template class SubsetEnumerator<StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>;
+        }
+    }
+}

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 <vector>
+
+
+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<typename DataType = std::nullptr_t>
+            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<uint_fast64_t> 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<uint_fast64_t> 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<k or the filter answers false in all cases).
+                bool setToFirstSubset();
+
+                //Increments the current subset.
+                //Returns true if there is a new subset and false if the current subset is already the last one.
+                bool incrementSubset();
+
+                // Default filter that returns always true.
+                static bool
+                trueFilter(std::vector<uint_fast64_t> 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<uint_fast64_t> current; //the current subset
+                std::vector<uint_fast64_t> 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_ */