Jonathan Kliem's profile - activity

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

The following answer should give you a step by step guide: https://ask.sagemath.org/question/354...

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

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

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
904076817914638457425434331054080000000

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

Sure. Thanks.

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")
True

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

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.

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")
True

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.

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

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()):
                it.ignore_subfaces()
        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):
                it.ignore_supfaces()
        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.

I have implemented this approach in https://trac.sagemath.org/ticket/29382.

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)
            data.pop(0)
            data.reverse()
            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:
                    continue
                for w in parse_indices(count, adjacencies[3:], self._cdd_H_to_sage_H):
                    if w is None:
                        continue
                    ret[v, w] = 1 
            ret.set_immutable()
            self.incidence_matrix.set_cache(ret)

        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)

Then

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

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}

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.