Accessing facets and affine hull data structures in C++

Questions and problems about using polymake go here.
Tryer
Posts: 13
Joined: 28 Feb 2022, 17:53

Accessing facets and affine hull data structures in C++

Postby Tryer » 03 Mar 2022, 03:41

Hello,

I have the following code in C++:

Code: Select all

#include <polymake/Main.h> #include <polymake/Matrix.h> #include <polymake/SparseMatrix.h> #include <polymake/Rational.h> using namespace polymake; int main(int argc, const char* argv[]) { try { const int dim = 4; Main pm; pm.set_application("polytope"); BigObject p("Polytope<Rational>"); std::vector<std::string> CoordinateLabels = { "x1", "x2", "x3" }; std::vector<std::vector<int>> Points; std::vector<int> point = { 0,0,0 }; Points.push_back(point); point = { 0,0,1 }; Points.push_back(point); point = { 0,1,0 }; Points.push_back(point); point = { 0,1,1 }; Points.push_back(point); point = { 1,0,0 }; Points.push_back(point); point = { 1,0,1 }; Points.push_back(point); point = { 1,1,0 }; Points.push_back(point); point = { 1,1,1 }; Points.push_back(point); p.take("VERTICES") << (ones_vector<Rational>() | 3*unit_matrix<Rational>(dim));//How to make p "take" or load the points from Points? const Matrix<Rational> f = p.give("FACETS"); cout << "facets" << endl << f << endl; const Matrix<Rational> affine_hull = p.give("AFFINE_HULL"); cout << "affine hull is " << endl << affine_hull<<endl; //How to access details of f and affine_hull object such as, how many facets are returned? In the example with Points, //is there something like f.size() which returns the number of facets? //f does not seem to have method size() implemented. Hence this query. //How should the user access the coefficient of x1 in the first facet, what is the rhs of the facet, //what is the sense of the facet ? //Essentially, can I iterate through the list of facets using something equivalent to an STL Iterator? //Similar query regarding accessing coefficients and right hand side of the affine hull equations. } catch (const std::exception& ex) { std::cerr << "ERROR: " << ex.what() << endl; return 1; } return 0; }
As specified in comments in the above code, I have the following queries.

(1)What is the syntax to populate the "VERTICES" property of a BigObject p("Polytope<Rational>") from the user's data structures such as say an array or std::vector of points? In the code example above, for instance, I have the eight vertices of a unit cube stored in a Points data structure that I would like to somehow pass into "VERTICES".

(2)Once polymake computes the facets and the affine hull, how can the user query details such as, how many facets / affine hull equations are there in the convex hull? f.size(), in the above code, for instance, does not compile as f does not seem to have the size() method. Similarly, how can one find out for the 3rd facet, what is the coefficient of variable x2, say, whether the facet inequality is a <= or a >= and what the rhs of the specific facet / affine hull is.

Thank you.

User avatar
gawrilow
Main Author
Posts: 405
Joined: 25 Dec 2010, 17:40

Re: Accessing facets and affine hull data structures in C++

Postby gawrilow » 03 Mar 2022, 09:34

BigObject operations give() and take() create objects which are conceptually similar to I/O streams. To set a property of a new object, you prepare it in an appropriate data structure and "stream" it into take(). Similarly, to retrieve a property, you "stream" the result of give() into an appropriate data structure, either using >> operator or copy initialization, like in your test program. To decide what can be the "appropriate" data structure for a given property, please refer to its documentation. Every property is defined along with a "canonical" type, which is used for internal storage and file serialization. For VERTICES, FACETS and AFFINE_HULL, this type is Matrix<Scalar>, where Scalar is the type parameter of the big object Polytope, that is, in your case, Rational. But you can also use structurally equivalent data structures, e.g. for a matrix it could be std::vector<std::vector<Scalar>>. Use of non-canonical data structures might be more convenient in some special cases, but comes with the price of hidden data conversion.

To prepare the matrix of VERTICES (or POINTS), fill it with coordinates, every row corresponds to one vertex. Beware of homogeneous coordinates! The concept is explained in tutorials dedicated to polyhedra, together with the encoding convention for facet inequalities.

Dimensions of a matrix can be queried as m.rows() resp. m.cols(), a single vector is retrieved with m.row(i), and a single element with m(i,j).

You might grep the source code in applications polytope or fan for give("VERTICES") or give("FACETS") for various "real life" examples of how matrices can be used; also check our wiki https://polymake.org/doku.php/user_guid ... nced_users .

User avatar
joswig
Main Author
Posts: 269
Joined: 24 Dec 2010, 11:10

Re: Accessing facets and affine hull data structures in C++

Postby joswig » 03 Mar 2022, 10:08

A few more remarks:

(1) Specifying polytopes in terms of a V-description works by providing POINTS. Using VERTICES is potentially more efficient, but this requires that the points (i.e., rows of the matrix) are actually, the vertices, and without repetitions.

(2) We use a homogeneous coordinate model. This also explains how inequalities are coded.

(3) You can browse polymake's source code on github to see how the give and take mechanism works. A good first example may be the function polarize which reads a polytope (or a cone) and writes its polar.

Note that the polymake client code looks a bit different than your code which uses libpolymake. However, this only affects "the outside" (i.e., headers, main, etc). That part of the code which actually does the work is the same.

Tryer
Posts: 13
Joined: 28 Feb 2022, 17:53

Re: Accessing facets and affine hull data structures in C++

Postby Tryer » 03 Mar 2022, 15:42

Thanks gawrilow and joswig. I have been able to make progress based on your comments. My code currently is:

Code: Select all

#include <polymake/Main.h> #include <polymake/Matrix.h> #include <polymake/SparseMatrix.h> #include <polymake/Rational.h> using namespace polymake; int main(int argc, const char* argv[]) { try { const int dim = 4; Main pm; pm.set_application("polytope"); BigObject p("Polytope<Rational>"); std::vector<std::string> CoordinateLabels = { "extravar", "x1", "x2", "x3" }; std::vector<std::vector<Int>> Points; std::vector<Int> point = { 1,0,0,0 }; Points.push_back(point); point = { 1,0,0,1 }; Points.push_back(point); point = { 1,0,1,0 }; Points.push_back(point); point = { 1,0,1,1 }; Points.push_back(point); point = { 1,1,0,0 }; Points.push_back(point); point = { 1,1,0,1 }; Points.push_back(point); point = { 1,1,1,0 }; Points.push_back(point); point = { 1,1,1,1 }; Points.push_back(point); p.take("VERTICES") << Points; const Matrix<Rational> f = p.give("FACETS"); cout << "facets" << endl << f << endl; cout << "Dimension of the facets matrix are " << f.rows() << " and " << f.cols() << endl; cout << "Printing the facets for easy reading" << endl; for (int i = 0; i < f.rows(); i++) { printf("Facet %d: ", i); Rational rhs = f(i, 0); for (int j = 1; j < f.cols(); j++) { Rational coef = f(i, j); if (coef.is_zero() == false) { cout << " + (" << coef << " " << CoordinateLabels[j] << ") "; } } cout << " >= - (" << rhs << ")" << endl; } const Matrix<Rational> affine_hull = p.give("AFFINE_HULL"); cout << "affine hull" << endl << affine_hull << endl; cout << "Dimension of the affine hull matrix are " << affine_hull.rows() << " and " << affine_hull.cols() << endl; cout << "Printing the affine hull for easy reading" << endl; for (int i = 0; i < affine_hull.rows(); i++) { printf("AffineHull %d: ", i); Rational rhs = affine_hull(i, 0); for (int j = 1; j < affine_hull.cols(); j++) { Rational coef = affine_hull(i, j); if (coef.is_zero() == false) { cout << " + (" << coef << " " << CoordinateLabels[j] << ") "; } } cout << " = - (" << rhs << ")" << endl; } } catch (const std::exception& ex) { std::cerr << "ERROR: " << ex.what() << endl; return 1; } return 0; }
and this produces output as :

Code: Select all

facets 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 -1 1 0 -1 0 1 -1 0 0 Dimension of the facets matrix are 6 and 4 Printing the facets for easy reading + (1 x3) >= - (0) + (1 x1) >= - (0) + (1 x2) >= - (0) + (-1 x3) >= - (1) + (-1 x2) >= - (1) + (-1 x1) >= - (1) affine hull Dimension of the affine hull matrix are 0 and 4 Printing the affine hull for easy reading
I have few specific questions:

(1)Is it safe to simply cast the Rational scalar to a double for further usage in my application?

That is, I currently have

Rational coef = affine_hull(i, j);

In my code, can I simply cast this as

double coef_dbl = (double)coef;

and use it subsequently without any loss of correctness?

(2)I currently incur the conversion hit that gawrilow alluded to since I use

std::vector<std::vector<Rational>> as opposed to Matrix<Rational>

as input for BigObject to "take" as points input.

What is the way to avoid this data conversion hit? Is there a Matrix<Rational> constructor that takes an std::vector<std::vector<Rational>> argument, or even better, std::vector<std::vector<double>> ?

I look at the documentation available at: https://polymake.org/doku.php/documenta ... mon#matrix

and it is not clear to me how to populate a Matrix<Rational> with my input data from my calling application.

ETA: I suppose that the best way here is to call clear(rows, cols) and then individually set each entry via function elem(r,c) since this returns an lvalue.

(3)Is Matrix<Rational> f 's access of element f(i,j) a constant time read / write operation?

User avatar
joswig
Main Author
Posts: 269
Joined: 24 Dec 2010, 11:10

Re: Accessing facets and affine hull data structures in C++

Postby joswig » 03 Mar 2022, 17:01

(1)Is it safe to simply cast the Rational scalar to a double for further usage in my application?
In general: no. Rational is an exact type, double is not. If a conversion makes sense or not will depend on what you want to compute.
Usually, you just want to stick with Rational. polymake supports arithmetic, linear algebra etc. Why do you want to convert at all?
(2) ETA: I suppose that the best way here is to call clear(rows, cols) and then individually set each entry via function elem(r,c) since this returns an lvalue.
Yes, this is one way.
For lots of examples check out https://github.com/polymake/polymake/bl ... ahedron.cc. This shows how to build matrices of type Matrix but also of the related types SparseMatrix and ListMatrix.
(3)Is Matrix<Rational> f 's access of element f(i,j) a constant time read / write operation?
Yes. For SparseMatrix access is logarithmic, for ListMatrix access is linear for the row index, and the cost for the columns depends on the choice of the row type. ListMatrix is useful if you want to build up a Matrix row by row, and you do not know how many rows you will have beforehand. Appending a new row is constant time (if the row already exists as a vector of the appropriate type).

Tryer
Posts: 13
Joined: 28 Feb 2022, 17:53

Re: Accessing facets and affine hull data structures in C++

Postby Tryer » 03 Mar 2022, 17:34

In general: no. Rational is an exact type, double is not. If a conversion makes sense or not will depend on what you want to compute.
Usually, you just want to stick with Rational. polymake supports arithmetic, linear algebra etc. Why do you want to convert at all?
I have other functions in my own code that expect double arguments. I have to use polymake's facet generation capabilities and then use this output to solve a different problem where all data is of double type. This part of my code has no knowledge of polymake or its data types. Hence the query about casting. Indeed, I believe I can make do with Matrix<float> and trade numerical precision for accuracy. All of my inputs to polymake will be rather small sized 0-1 vectors (around 6 dimensions) with some special structure.

User avatar
joswig
Main Author
Posts: 269
Joined: 24 Dec 2010, 11:10

Re: Accessing facets and affine hull data structures in C++

Postby joswig » 03 Mar 2022, 18:42

There are several things to consider then:

(1) If you pass float coordinates as point coordinates of a Polytope object, then they will be converted to exact rationals automatically.

(2) In general, there is no way (known) to "approximately" compute convex hulls with floats throughout. In dimension 3 one can do this (with a bit of effort), but not in dimensions four or higher. While there exists code (like, e.g., qhull) which does compute convex hulls with floats, it may happen that the (exact or approximate) incidences of vertices and facets do not fit any polytope.

(3) Stick with exact computations as long as possible. Convert at the very last moment.

(4) If you are computing in six dimensions only it may not matter that much, but it is worth knowing that various convex hull codes behave very differently with respect to speed on various input; see https://link.springer.com/article/10.10 ... 016-0104-z. For 0/1-polytopes polymake's default (which is ppl) should be good.


Return to “Helpdesk”