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/45538/display-vertex-label-a-polytope/ 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 included?