Labelling vertices in a displayed polytope

edit

I am working with the following code:

def make_polytope(v,d):
    n = len(v)

    I = identity_matrix(n)
    zeros_column = matrix(ZZ, n, 1, [0]*n)
    matrix_with_zeros = zeros_column.augment(I)

    vector1 = vector([-d] + v.list())
    vector2 = vector([d] + (-v).list())

    M = matrix_with_zeros.stack(matrix([vector1, vector2]))

    P1 = Polyhedron(ieqs = M)

    P = Polyhedron(vertices = P1.integral_points())

    return P
a = vector([5,19,27,31])
d = 81
P = make_polytope(a,d)
i=0
for f in P.vertices():
    print(i,": ", f)
    i=i+1
plot = P.plot(line='red', polygon=False, point={'size':30,'color':'green'}, vertex_labels=True)

plot.show()

If you run it, you get a polytope in 3D space printed out (changing a and d gives you different polytopes, but this is both simple and complicated enough for this purpose).

Note that with the above data, P is a 3D polytope in 4-space. As I understand it, Sage notices this and projects P to 3D space "orthonormally" via plot(**) and then prints it out with show().

It's all very impressive but I am not able to better study the polytope, as it does not come with labels. Whether I keep or remove the vertex_labels command, it makes no difference.

I have read both this: https://ask.sagemath.org/question/455... and the comment and link within and it is not helpful (neither it seems to provide a solution, nor I want this done in tikZ for the moment).

Q: How can I change the above for labels for vertices to be plotted?