19template<
typename ValueType>
21 if (halfspaces.empty()) {
23 emptyStatus = EmptyStatus::Nonempty;
25 Eigen::Index maxCol = halfspaces.front().normalVector().size();
26 Eigen::Index maxRow = halfspaces.size();
29 for (Eigen::Index row = 0; row < A.rows(); ++row) {
30 assert((Eigen::Index)halfspaces[row].normalVector().size() == maxCol);
31 b(row) = halfspaces[row].offset();
34 emptyStatus = EmptyStatus::Unknown;
38template<
typename ValueType>
41 emptyStatus = EmptyStatus::Empty;
43 std::vector<EigenVector> eigenPoints;
44 eigenPoints.reserve(points.size());
45 for (
auto const& p : points) {
53 emptyStatus = EmptyStatus::Nonempty;
57template<
typename ValueType>
59 boost::optional<std::vector<Point>>
const& points) {
61 STORM_LOG_WARN_COND(!points,
"Creating a NativePolytope where halfspaces AND points are given. The points will be ignored.");
62 return std::make_shared<NativePolytope<ValueType>>(*halfspaces);
64 return std::make_shared<NativePolytope<ValueType>>(*points);
66 STORM_LOG_THROW(
false, storm::exceptions::UnexpectedException,
"Creating a NativePolytope but no representation was given.");
70template<
typename ValueType>
75template<
typename ValueType>
77 : emptyStatus(
std::move(other.emptyStatus)), A(
std::move(other.A)), b(
std::move(other.b)) {
81template<
typename ValueType>
83 : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) {
87template<
typename ValueType>
89 : emptyStatus(emptyStatus), A(halfspaceMatrix), b(halfspaceVector) {
93template<
typename ValueType>
98template<
typename ValueType>
100 std::vector<EigenVector> eigenVertices = getEigenVertices();
101 std::vector<Point> result;
102 result.reserve(eigenVertices.size());
103 for (
auto const& p : eigenVertices) {
109template<
typename ValueType>
111 std::vector<Halfspace<ValueType>> result;
112 result.reserve(A.rows());
114 for (Eigen::Index row = 0; row < A.rows(); ++row) {
120template<
typename ValueType>
122 if (emptyStatus == EmptyStatus::Unknown) {
125 std::vector<storm::expressions::Expression> constraints = getConstraints(*manager, declareVariables(*manager,
"x"));
126 for (
auto const& constraint : constraints) {
127 solver->add(constraint);
129 switch (solver->check()) {
131 emptyStatus = EmptyStatus::Nonempty;
134 emptyStatus = EmptyStatus::Empty;
137 STORM_LOG_THROW(
false, storm::exceptions::UnexpectedException,
"Unexpected result of SMT solver during emptyness-check of Polytope.");
141 return emptyStatus == EmptyStatus::Empty;
144template<
typename ValueType>
146 return A.rows() == 0;
149template<
typename ValueType>
152 for (Eigen::Index row = 0; row < A.rows(); ++row) {
153 if ((A.row(row) * x)(0) > b(row)) {
160template<
typename ValueType>
162 STORM_LOG_THROW(other->isNativePolytope(), storm::exceptions::InvalidArgumentException,
163 "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported");
164 if (this->isUniversal()) {
166 }
else if (other->isUniversal()) {
172 std::vector<storm::expressions::Variable> variables = declareVariables(*manager,
"x");
173 std::vector<storm::expressions::Expression> constraints = getConstraints(*manager, variables);
175 for (
auto const& constraint : constraints) {
176 constraintsThis = constraintsThis && constraint;
178 solver->add(!constraintsThis);
180 for (
auto const& constraint : constraints) {
181 solver->add(constraint);
183 switch (solver->check()) {
189 STORM_LOG_THROW(
false, storm::exceptions::UnexpectedException,
"Unexpected result of SMT solver during containment check of two polytopes.");
195template<
typename ValueType>
197 STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException,
198 "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported");
199 if (this->isUniversal()) {
201 }
else if (rhs->isUniversal()) {
202 return std::make_shared<NativePolytope<ValueType>>(*this);
205 EigenMatrix resultA(A.rows() + nativeRhs.A.rows(), A.cols());
206 resultA << A, nativeRhs.A;
208 resultb << b, nativeRhs.b;
209 return std::make_shared<NativePolytope<ValueType>>(EmptyStatus::Unknown, std::move(resultA), std::move(resultb));
213template<
typename ValueType>
219 resultb(0) = halfspace.
offset();
220 return std::make_shared<NativePolytope<ValueType>>(EmptyStatus::Unknown, std::move(resultA), std::move(resultb));
225 resultb << b, halfspace.
offset();
226 return std::make_shared<NativePolytope<ValueType>>(EmptyStatus::Unknown, std::move(resultA), std::move(resultb));
229template<
typename ValueType>
231 STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException,
232 "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported");
233 if (this->isEmpty()) {
235 }
else if (rhs->isEmpty()) {
236 return std::make_shared<NativePolytope<ValueType>>(*this);
237 }
else if (this->isUniversal() || rhs->isUniversal()) {
238 return std::make_shared<NativePolytope<ValueType>>(std::vector<Halfspace<ValueType>>());
241 STORM_LOG_WARN_COND_DEBUG(
false,
"Implementation of convex union of two polytopes only works if the polytopes are bounded. This is not checked.");
244 std::vector<EigenVector> resultVertices = this->getEigenVertices();
245 resultVertices.insert(resultVertices.end(), std::make_move_iterator(rhsVertices.begin()), std::make_move_iterator(rhsVertices.end()));
252template<
typename ValueType>
254 STORM_LOG_THROW(rhs->isNativePolytope(), storm::exceptions::InvalidArgumentException,
255 "Invoked operation between a NativePolytope and a different polytope implementation. This is not supported");
258 if (this->isEmpty() || nativeRhs.
isEmpty()) {
259 return std::make_shared<NativePolytope<ValueType>>(std::vector<Point>());
262 std::vector<std::pair<EigenVector, ValueType>> resultConstraints;
263 resultConstraints.reserve(A.rows() + nativeRhs.A.rows());
266 for (Eigen::Index
i = 0;
i < A.rows(); ++
i) {
267 auto optimizationRes = nativeRhs.
optimize(A.row(
i));
268 if (optimizationRes.second) {
269 resultConstraints.emplace_back(A.row(
i), b(
i) + (A.row(
i) * optimizationRes.first)(0));
275 for (Eigen::Index
i = 0;
i < nativeRhs.A.rows(); ++
i) {
276 auto optimizationRes = optimize(nativeRhs.A.row(
i));
277 if (optimizationRes.second) {
278 resultConstraints.emplace_back(nativeRhs.A.row(
i), nativeRhs.b(
i) + (nativeRhs.A.row(
i) * optimizationRes.first)(0));
283 if (resultConstraints.empty()) {
284 return std::make_shared<NativePolytope<ValueType>>(std::vector<Halfspace<ValueType>>());
286 EigenMatrix newA(resultConstraints.size(), resultConstraints.front().first.rows());
288 for (Eigen::Index
i = 0;
i < newA.rows(); ++
i) {
289 newA.row(
i) = resultConstraints[
i].first;
290 newb(
i) = resultConstraints[
i].second;
292 return std::make_shared<NativePolytope<ValueType>>(EmptyStatus::Nonempty, std::move(newA), std::move(newb));
296template<
typename ValueType>
298 STORM_LOG_THROW(!matrix.empty(), storm::exceptions::InvalidArgumentException,
"Invoked affine transformation with a matrix without rows.");
299 Eigen::Index rows = matrix.size();
300 Eigen::Index columns = matrix.front().size();
302 for (Eigen::Index row = 0; row < rows; ++row) {
307 Eigen::FullPivLU<EigenMatrix> luMatrix(eigenMatrix);
308 STORM_LOG_THROW(luMatrix.isInvertible(), storm::exceptions::NotImplementedException,
309 "Affine Transformation of native polytope only implemented if the transformation matrix is invertable");
311 return std::make_shared<NativePolytope<ValueType>>(std::vector<Halfspace<ValueType>>());
315 return std::make_shared<NativePolytope<ValueType>>(emptyStatus, std::move(newA), std::move(newb));
318template<
typename ValueType>
321 return std::make_pair(
Point(),
false);
325 std::vector<storm::expressions::Variable> variables;
326 variables.reserve(A.cols());
327 for (Eigen::Index
i = 0;
i < A.cols(); ++
i) {
330 std::vector<storm::expressions::Expression> constraints = getConstraints(solver.
getManager(), variables);
331 for (
auto const& constraint : constraints) {
337 auto result = std::make_pair(
Point(),
true);
338 result.first.reserve(variables.size());
339 for (
auto const& var : variables) {
345 return std::make_pair(
Point(),
false);
349template<
typename ValueType>
352 return std::make_pair(EigenVector(),
false);
356 std::vector<storm::expressions::Variable> variables;
357 variables.reserve(A.cols());
358 for (Eigen::Index
i = 0;
i < A.cols(); ++
i) {
359 variables.push_back(solver.addUnboundedContinuousVariable(
"x" + std::to_string(
i),
static_cast<ValueType
>(direction(
i))));
361 std::vector<storm::expressions::Expression> constraints = getConstraints(solver.getManager(), variables);
362 for (
auto const& constraint : constraints) {
363 solver.addConstraint(
"", constraint);
367 if (solver.isOptimal()) {
368 auto result = std::make_pair(EigenVector(A.cols()),
true);
369 for (Eigen::Index
i = 0;
i < A.cols(); ++
i) {
370 result.first(
i) = solver.getContinuousValue(variables[
i]);
375 return std::make_pair(EigenVector(),
false);
379template<
typename ValueType>
383template<
typename ValueType>
390template<
typename ValueType>
392 std::string
const& namePrefix)
const {
393 std::vector<storm::expressions::Variable> result;
394 result.reserve(A.cols());
395 for (Eigen::Index col = 0; col < A.cols(); ++col) {
396 result.push_back(manager.declareVariable(namePrefix + std::to_string(col), manager.getRationalType()));
401template<
typename ValueType>
403 std::vector<storm::expressions::Variable>
const& variables)
const {
404 std::vector<storm::expressions::Expression> result;
405 for (Eigen::Index row = 0; row < A.rows(); ++row) {
407 for (Eigen::Index col = 1; col < A.cols(); ++col) {
408 lhs = lhs + manager.rational(A(row, col)) * variables[col].getExpression();
410 result.push_back(lhs <= manager.rational(b(row)));
415template<
typename ValueType>
418 return create(boost::none, {});
422 auto variables = declareVariables(*manager,
"x");
423 std::vector<storm::expressions::Expression> constraints = getConstraints(*manager, variables);
424 for (
auto const& constraint : constraints) {
425 solver->add(constraint);
428 for (Eigen::Index row = 0; row < A.rows(); ++row) {
430 for (Eigen::Index col = 1; col < A.cols(); ++col) {
431 lhs = lhs + manager->rational(A(row, col)) * variables[col].getExpression();
434 solver->add(lhs >= manager->rational(b(row)));
435 switch (solver->check()) {
437 keptConstraints.
set(row,
true);
441 STORM_LOG_THROW(
false, storm::exceptions::UnexpectedException,
"Unexpected result of SMT solver during emptyness-check of Polytope.");
446 std::vector<Halfspace<ValueType>> newHalfspaces;
448 for (
auto row : keptConstraints) {
451 return create(newHalfspaces, boost::none);
455#ifdef STORM_HAVE_CARL
static Eigen::Matrix< ValueType, Eigen::Dynamic, 1 > toEigenVector(std::vector< ValueType > const &v)
static std::vector< ValueType > toStdVector(Eigen::Matrix< ValueType, Eigen::Dynamic, 1 > const &v)
This class is responsible for managing a set of typed variables and all expressions using these varia...
storm::expressions::ExpressionManager const & getManager() const
Retrieves the manager for the variables created for this solver.
Variable addUnboundedContinuousVariable(std::string const &name, ValueType objectiveFunctionCoefficient=0)
Registers a unbounded continuous variable, i.e.
A class that implements the LpSolver interface using Z3.
virtual void update() const override
Updates the model to make the variables that have been declared since the last call to update usable.
virtual ValueType getContinuousValue(Variable const &variable) const override
Retrieves the value of the continuous variable with the given name.
virtual void addConstraint(std::string const &name, Constraint const &constraint) override
Adds a the given constraint to the LP problem.
virtual bool isOptimal() const override
Retrieves whether the model was found to be optimal, i.e.
virtual void optimize() const override
Optimizes the LP problem previously constructed.
A bit vector that is internally represented as a vector of 64-bit values.
void set(uint_fast64_t index, bool value=true)
Sets the given truth value at the given index.
uint_fast64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
std::vector< ValueType > const & normalVector() const
ValueType const & offset() const
std::vector< EigenVector > & getResultVertices()
void generateVerticesFromConstraints(EigenMatrix const &constraintMatrix, EigenVector const &constraintVector, bool generateRelevantHyperplanesAndVertexSets)
virtual bool contains(Point const &point) const override
Returns true iff the given point is inside of the polytope.
NativePolytope(std::vector< Halfspace< ValueType > > const &halfspaces)
Creates a NativePolytope from the given halfspaces The resulting polytope is defined as the intersect...
Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic > EigenMatrix
virtual std::vector< storm::expressions::Variable > declareVariables(storm::expressions::ExpressionManager &manager, std::string const &namePrefix) const override
declares one variable for each dimension and returns the obtained variables.
virtual std::pair< Point, bool > optimize(Point const &direction) const override
Finds an optimal point inside this polytope w.r.t.
virtual std::shared_ptr< Polytope< ValueType > > clean() override
Performs cleaning operations, e.g., deleting redundant halfspaces.
virtual std::shared_ptr< Polytope< ValueType > > affineTransformation(std::vector< Point > const &matrix, Point const &vector) const override
Returns the affine transformation of this polytope P w.r.t.
virtual bool isNativePolytope() const override
Polytope< ValueType >::Point Point
virtual std::vector< Point > getVertices() const override
Returns the vertices of this polytope.
Eigen::Matrix< ValueType, Eigen::Dynamic, 1 > EigenVector
virtual bool isUniversal() const override
Returns whether this polytope is universal (i.e., equals R^n).
virtual std::shared_ptr< Polytope< ValueType > > intersection(std::shared_ptr< Polytope< ValueType > > const &rhs) const override
Intersects this polytope with rhs and returns the result.
virtual std::shared_ptr< Polytope< ValueType > > minkowskiSum(std::shared_ptr< Polytope< ValueType > > const &rhs) const override
Returns the minkowskiSum of this polytope and rhs.
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 or points.
virtual bool isEmpty() const override
Returns whether this polytope is the empty set.
virtual std::vector< storm::expressions::Expression > getConstraints(storm::expressions::ExpressionManager const &manager, std::vector< storm::expressions::Variable > const &variables) const override
returns the constrains defined by this polytope as an expression over the given variables
virtual ~NativePolytope()
virtual std::vector< Halfspace< ValueType > > getHalfspaces() const override
Returns the halfspaces of this polytope.
virtual std::shared_ptr< Polytope< ValueType > > convexUnion(std::shared_ptr< Polytope< ValueType > > const &rhs) const override
Returns the convex union of this polytope and rhs.
void generateHalfspacesFromPoints(std::vector< EigenVector > &points, bool generateRelevantVerticesAndVertexSets)
EigenMatrix & getResultMatrix()
EigenVector & getResultVector()
virtual std::unique_ptr< storm::solver::SmtSolver > create(storm::expressions::ExpressionManager &manager) const
Creates a new SMT solver instance.
#define STORM_LOG_WARN_COND(cond, message)
#define STORM_LOG_THROW(cond, exception, message)
#define STORM_LOG_WARN_COND_DEBUG(cond, message)