Ask Your Question

Revision history [back]

As far as i know, there is no canonical measure on the set of polytopes with a given number of vertices. So, you need to provide what "random" means for you. Otherwise, here is a basic rejection algorithm:

def random_polytope(npoints=4, dimension=3, ring=RDF):
    while True:
        P = Polyhedron(vertices=[(ring.random_element() for i in range(dimension)) for n in range(npoints)])
        if P.n_vertices() == npoints:
            return P

Then, you can do:

sage: random_polytope(10,3)
A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 10 vertices

Note that when the dimension is small and the number of points is big, you might wait a long time since in most cases, one of the points will be contained in the convex hull of the others.

As far as i know, there is no canonical measure on the set of polytopes with a given number of vertices. So, you need to provide what "random" means for you. Otherwise, here is a basic rejection algorithm:

def random_polytope(npoints=4, dimension=3, ring=RDF):
    while True:
        P = Polyhedron(vertices=[(ring.random_element() for i in range(dimension)) for n in range(npoints)])
        if P.n_vertices() == npoints:
            return P

Then, you can do:

sage: random_polytope(10,3)
A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 10 vertices

Note that when the dimension is small and the number of points is big, you might wait a long time since in most cases, one of the points will be contained in the convex hull of the others.

EDIT

Here is how to make a polytope from points uniformly sampled on the unit sphere: it suffice to take vectors whose coordinates are sampled independently accorfing to the normal distribution and then project them on the sphere:

def random_polytope(npoints=4, dimension=3):
    vert = []
    for i in range(npoints):
        v = vector((gauss(0,1) for d in range(dimension)))
        v = v / v.norm()
        vert.append(v)
    return Polyhedron(vertices=vert)

Then you can do:

sage: p = random_polytope(npoints=100, dimension=3)
sage: p.plot()
Launched jmol viewer for Graphics3d Object

If you want to sample simple d-polytope from this, you can use a rejection algorithm as before.