Ask Your Question

Jonathan Kliem's profile - activity

2021-03-25 19:50:49 +0200 commented question Oriented vertices of Polyhedron

Am I understanding correctly that you are referring only to 3-dimensional polytopes? (Because I woudn't even understand

2020-08-25 10:59:18 +0200 commented answer Computing projection of a polyhedron

The following answer should give you a step by step guide:

In particular how to obtain the matrix that induces a projection to an affine span.

2020-08-24 08:18:30 +0200 commented answer Computing projection of a polyhedron

In sage 9.1 one can skip the third part by just multiplying the matrix with the polyhedron: proj_mat*P.

2020-08-13 18:46:13 +0200 received badge  Necromancer (source)
2020-07-15 18:18:51 +0200 answered a question Parallelizing Polyhedral volume computation

In sage 9.2+ (and probably since a while), you can use the optional package pynormaliz to do parallel volume computations. If you have sage compiled with

MAKE='make -j64'

one a computer with 64 processors you should use that many processors (I hope at least, it worked for me for 8).

The following used 8 threads for me:

sage: P = polytopes.cyclic_polytope(10,20, backend='normaliz')   # optional - pynormaliz
sage: P.volume()  # optional pynormaliz
2020-05-11 21:54:17 +0200 received badge  Nice Answer (source)
2020-05-11 21:16:32 +0200 commented answer What is the most efficient way to "look up" a face in the face lattice of a polyhedron?

I edited the answer to explain what the code is doing.

2020-04-24 09:44:34 +0200 commented answer Projection along affine hull

Sure. Thanks.

2020-04-24 09:41:12 +0200 received badge  Necromancer (source)
2020-04-23 23:30:25 +0200 commented question Can sage determine if a cone is gorenstein?

Normaliz might be useful here. It's an optional package of sage and if you have it installed you can use it:

sage: import PyNormaliz
sage: C = PyNormaliz.NmzCone("cone", [[0, 0, 1, -1], [0, 1, -1, 0], [1, -1, 0, 0]])  # homogenous input
sage: PyNormaliz.NmzResult(C, "IsGorenstein")

I'm not sure, whether or not the definition in normaliz agrees with your definition.

2020-04-23 23:05:13 +0200 received badge  Necromancer (source)
2020-04-23 22:09:03 +0200 commented answer What is the most efficient way to "look up" a face in the face lattice of a polyhedron?

Yes, I noticed. It wasn't until JP pointed me to this post that I realized that this would be a nice thing to have. This is all very much work in progress and it reminds me that more stuff of combinatorial polyhedron should be exposed. I might be about the only person aware of the methods ignore_subfaces and ignore_supfaces.

2020-04-23 21:43:23 +0200 answered a question Can sage determine if a cone is gorenstein?

Normaliz can do it, which is an optional package of sage. So I guess the following answer is only useful, if you have installed sage from source. (Technically, all you need is a normaliz and a way to transfer your input into normaliz input.)

Unfortunately, normaliz only accepts homogeneous format for this property and as of now (sage 9.1), we have to do some work :

sage: import PyNormaliz
sage: P = Polyhedron(rays=[[0, 0, 1, -1], [0, 1, -1, 0], [1/2, -1/2, 0, 0]], backend='normaliz', verbose=True)
# ----8<---- Equivalent Normaliz input file ----8<---- 
amb_space 4
cone 3
 0 0 1 -1
 0 1 -1 0
 1 -1 0 0
subspace 0
vertices 1
 0 0 0 0 1
# ----8<-------------------8<-------------------8<----
# Calling PyNormaliz.NmzCone(cone=[[0, 0, 1, -1], [0, 1, -1, 0], [1, -1, 0, 0]], subspace=[], vertices=[[0, 0, 0, 0, 1]])
sage: C = PyNormaliz.NmzCone("cone", [[0, 0, 1, -1, 0], [0, 1, -1, 0, 0], [1, -1, 0, 0, 0], [0, 0, 0, 0, 1]])  # homogenized input
sage: PyNormaliz.NmzResult(C, "IsGorenstein")

The verbose call of the polyhedron constructor is of course not necessary, but it is useful to figure out the normaliz input. The homogenized input is obtained by appending a 0 to rays and a 1 to vertices and putting them together as input cone.

2020-04-23 18:54:21 +0200 answered a question Projection along affine hull

In sage 9.1+ matrices act on polyhedra as linear transformation. This makes it easy.

Let me explain how to project the 3-dimensional permutahedron P to the affine span of [(3, 2, 4, 1), (4, 1, 3, 2), (1, 2, 3, 4)]. By translation, we can reduce the problem to a linear one. We can instead linearly project P1 = P - vector((3, 2, 4, 1)) to the linear span of [((1, -1, -1, 1), (-2, 0, -1, 3))].

Obtaining P1:

sage: P = polytopes.permutahedron(4)  # dim 3 and ambient dim 4
sage: P1 = P - vector((3, 2, 4, 1))

To do this we use the right kernel:

sage: M = Matrix(((1, -1, -1, 1), (-2, 0, -1, 3)))
sage: B = M.stack(M.right_kernel_matrix()).transpose()
sage: B
[ 1 -2  1  0]
[-1  0  1  2]
[-1 -1  1 -3]
[ 1  3  1 -1]

The first two rows are spanning our linear space we want to map to. The third and fourth row are a base extension such that they are orthogonal to the first and second row.

Now we want to linearly transform P such that the vectors corresponding to the first and second row get mapped to themselves and the others to zero. That's linear algebra:

sage: P2 = B*Matrix([[1,0,0,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])*B.inverse()*P1

Finally, we add the vector again:

sage: P2 + vector((3, 2, 4, 1))
A 2-dimensional polyhedron in QQ^4 defined as the convex hull of 10 vertices
2020-04-23 18:04:54 +0200 answered a question What is the most efficient way to "look up" a face in the face lattice of a polyhedron?

This is a very natural question, but unfortunately at the moment the answer is a bit complicated. This answer assumes sage 9.1+:

As this is a combinatorial property, you can go into CombinatorialPolyhedron, which is now a method of a polyhedron. There you can obtain the answer via the face iterator. The face iterator iterates through all faces of the polyhedron with a depth-first search. Of course you could just iterate until you what you are looking for, but the iterator can also be manipulated. For this I need to explain a bit, how it works. Assume for now that the iterator is in non-dual mode, it generates all faces from the facets.

The iterator yields first all facets. Then it inductively yields all faces of the first facet (treating it as it's own polyhedron) and then marks it as being visited. Then it yields all faces of the second facet, but the faces contained in the first facet. After that it marks the second facet as being visited and so on.

Now in our situation it might happen, that the first facet does not contain the goal face. Then we want to mark it immediately as being visited. Likewise we proceed with all facets. The below code will proceed like this for the ridges and so on, which might be easier to adapt, when you need different properties.

You might observe that in this particular situation it suffices to mark the correct facets as visited.

Now, if your polyhedron contains many facets but few vertices, generating all faces from the facets might be a bad idea. This is why the iterator can work in dual mode as well. In this case you can decide to skip supfaces instead of subfaces.

As mentioned before, this is a natural question and eventually I think there should be a method like find-join-of-vertices or find-meet-of-facets.

The following should work (it assumes that the face really does exist, otherwise you need to alter the termination condition of the while loop to fit your need):

from sage.geometry.polyhedron.face import combinatorial_face_to_polyhedral_face
def is_subset(a,b):
    return all(x in b for x in a)

def obtain_face(P, *indices):
    indices = tuple(indices)
    C = P.combinatorial_polyhedron()
    it = C.face_iter()

    # Let f be the first facet in non-dual mode or the first vertex in dual mode.
    f = next(it)

    while f.ambient_V_indices() != indices:
        if not it.dual:
            # Skipping faces unless our goal face is a subface.
            if not is_subset(indices, f.ambient_V_indices()):
        else: # the face iterator starts with the vertices
            # Skipping faces unless our goal face is a supface.
            if not is_subset(f.ambient_V_indices(), indices):
        f = next(it)

    # f is now the combinatorial face we where looking for (assuming it exists).
    # Obtain a polyhedral face from it.
    return combinatorial_face_to_polyhedral_face(P, f)

Then you get:

sage: P = polytopes.cross_polytope(4)
sage: f = obtain_face(P, 2,4,6)
sage: f.ambient_V_indices()
(2, 4, 6)
sage: f
A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 3 vertices
sage: P = polytopes.hypercube(3)
sage: f = obtain_face(P, 0,1,5,6)
sage: f.ambient_V_indices()
(0, 1, 5, 6)
sage: f
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices

This should be pretty fast. Note, that no options of the face iterator are exposed in the polyhedron and one has to go into combinatorial polyhedron (skipping subfaces or supfaces is crucial, because we only want one depth search and find our face immediately).

The the subset check function is not optimized deliberately.

2020-03-27 11:45:05 +0200 commented answer Compute the centroid of a polytope

I have implemented this approach in

2020-03-20 12:31:42 +0200 received badge  Nice Answer (source)
2020-03-19 18:45:10 +0200 received badge  Necromancer (source)
2020-03-19 16:46:15 +0200 answered a question VoronoiDiagram returns empty regions

I think I have fixed the bug in ticket 29176.

Hopefully, it will be available from Sage version 9.1+.

In the meantime you can use the following workaround.

from sage.all import *
from sage.geometry.polyhedron.parent import Polyhedra_RDF_cdd
from sage.geometry.polyhedron.backend_cdd import Polyhedron_RDF_cdd
class Polyhedron_RDF_cdd_workaround(Polyhedron_RDF_cdd):
    def _init_from_cdd_output(self, cddout):
        Polyhedron_RDF_cdd._init_from_cdd_output(self, cddout)
        cddout = cddout.splitlines()

        def parse_indices(count, cdd_indices, cdd_indices_to_sage_indices=None):
            cdd_indices = [int(x) for x in cdd_indices]
            if cdd_indices_to_sage_indices is None:
                cdd_indices_to_sage_indices = {i: i-1 for i in cdd_indices}
            if count < 0:
                assert cdd_indices_to_sage_indices is not None, "Did not expect negative counts here"
                count = -count
                cdd_indices = list(set(cdd_indices_to_sage_indices.keys()) - set(cdd_indices))
                assert count in [len(cdd_indices), len(cdd_indices) - 1]
            assert count == len(cdd_indices)
            return [cdd_indices_to_sage_indices[i] for i in cdd_indices if cdd_indices_to_sage_indices[i] is not None]

        def parse_incidence_matrix(intro, data):
            if 'incidence_matrix' in self.__dict__:
                raise NotImplementedError("can not replace internal representation as this breaks caching")
            N = len(self._Hrepresentation)
            M = len(self._Vrepresentation)
            ret = matrix(ZZ, M, N, 0)
            for adjacencies in data:
                assert adjacencies[2] == ':', "Not a line of adjacency data"
                cdd_vertex = int(adjacencies[0])
                count = int(adjacencies[1])
                if cdd_vertex not in self._cdd_V_to_sage_V:
                    self._cdd_V_to_sage_V[cdd_vertex] = None
                v = self._cdd_V_to_sage_V[cdd_vertex]
                if v is None:
                for w in parse_indices(count, adjacencies[3:], self._cdd_H_to_sage_H):
                    if w is None:
                    ret[v, w] = 1 

        Polyhedron_RDF_cdd._parse_block(cddout, 'Vertex incidence', parse_incidence_matrix)

class Polyhedra_RDF_cdd_workaround(Polyhedra_RDF_cdd):
    Element = Polyhedron_RDF_cdd_workaround

class VoronoiDiagram_workaround(VoronoiDiagram):
    def __init__(self, points):
        VoronoiDiagram.__init__(self, points)
        self._P = {}
        if self._n > 0:
            e = [([sum(vector(i)[k] ** 2
                   for k in range(self._d))] +
                  [(-2) * vector(i)[l] for l in range(self._d)] + [1])
                 for i in self._points]
            e = [[self._base_ring(i) for i in k] for k in e]
            parent = Polyhedra_RDF_cdd_workaround(RDF, self._d + 1, 'cdd')
            p = parent(None, [e, []])
        for i in range(self._n):
            parent = Polyhedra_RDF_cdd_workaround(RDF, self._d, 'cdd')
            equ = p.Hrepresentation(i)
            pvert = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_vertex()]
            prays = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_ray()]
            pline = [[u[k] for k in range(self._d)] for u in equ.incident()
                     if u.is_line()]
            (self._P)[self._points[i]] = parent([pvert, prays, pline], None)


P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
V = VoronoiDiagram_workaround(P)
R = V.regions()

will produce correctly

{P(-2687.19000000000, -2088.53000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line,
 P(-2686.81000000000, -2084.19000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line}
2019-08-27 22:29:02 +0200 received badge  Supporter (source)
2019-08-27 22:24:33 +0200 received badge  Editor (source)
2019-08-27 22:11:31 +0200 received badge  Necromancer (source)
2019-08-27 22:11:31 +0200 received badge  Teacher (source)
2019-08-27 19:08:27 +0200 answered a question Forming a polytope from only its combinatorial data

From Sage version 8.9+ there will be CombinatorialPolyhedron, which can be initialized with a list of facets, as

CombinatorialPolyhedron([[1,2], [2,3], [3,4], [4,1]])

With this object you have a face iterator, the face lattice and some other methods. I don't know if that helps.