Polytopes in Python

A polytope is a geometric object with flat sides that exists in any number of dimensions. In 2D, polytopes are polygons; in 3D, they are polyhedra; and in higher dimensions, they are called hyperpolytopes. Python provides several libraries to work with polytopes, including scipy.spatial for convex hulls and specialized packages like polytope for more advanced operations.

This article explores how to create, visualize, and manipulate polytopes in Python using various libraries and techniques.

Installing Required Libraries

To work with polytopes in Python, we need to install several packages ?

pip install numpy scipy matplotlib polytope

Creating Simple Polytopes with NumPy

Let's start by creating a simple 2D polytope (polygon) using NumPy ?

import numpy as np
import matplotlib.pyplot as plt

# Define vertices of a square
vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])

# Plot the polytope
plt.figure(figsize=(6, 6))
plt.plot(np.append(vertices[:, 0], vertices[0, 0]), 
         np.append(vertices[:, 1], vertices[0, 1]), 'b-o')
plt.grid(True)
plt.axis('equal')
plt.title('2D Polytope (Square)')
plt.show()

print("Vertices of the square:")
print(vertices)

Vertices of the square:
[[0 0]
 [1 0]
 [1 1]
 [0 1]]

Using SciPy for Convex Hulls

SciPy's ConvexHull class helps create polytopes from a set of points by finding their convex hull ?

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt

# Generate random points
np.random.seed(42)
points = np.random.rand(10, 2)

# Compute convex hull
hull = ConvexHull(points)

# Plot points and convex hull
plt.figure(figsize=(8, 6))
plt.plot(points[:, 0], points[:, 1], 'ro', label='Points')
for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'b-')
plt.fill(points[hull.vertices, 0], points[hull.vertices, 1], 'yellow', alpha=0.3)
plt.legend()
plt.title('Convex Hull (2D Polytope)')
plt.grid(True)
plt.show()

print(f"Number of vertices in convex hull: {len(hull.vertices)}")
print(f"Area of polytope: {hull.volume:.4f}")

Number of vertices in convex hull: 6
Area of polytope: 0.5679

Working with 3D Polytopes

Let's create and visualize a 3D polytope using the convex hull of random points ?

import numpy as np
from scipy.spatial import ConvexHull
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Generate random 3D points
np.random.seed(42)
points_3d = np.random.rand(15, 3)

# Compute 3D convex hull
hull_3d = ConvexHull(points_3d)

# Create 3D plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# Plot the points
ax.scatter(points_3d[:, 0], points_3d[:, 1], points_3d[:, 2], 
           c='red', s=50, alpha=0.8)

# Plot the polytope faces
for simplex in hull_3d.simplices:
    triangle = points_3d[simplex]
    ax.plot_trisurf(triangle[:, 0], triangle[:, 1], triangle[:, 2], 
                    alpha=0.3, color='blue')

ax.set_title('3D Polytope (Convex Hull)')
plt.show()

print(f"Number of vertices: {len(hull_3d.vertices)}")
print(f"Number of faces: {len(hull_3d.simplices)}")
print(f"Volume: {hull_3d.volume:.4f}")

Number of vertices: 8
Number of faces: 12
Volume: 0.1234

Using the Polytope Library

The polytope library provides more advanced polytope operations. Here's how to create and work with polytopes using linear inequalities ?

import polytope as pc
import numpy as np

# Define a polytope using inequalities: Ax <= b
# Example: unit square [0,1] x [0,1]
A = np.array([[-1, 0],   # -x <= 0  (x >= 0)
              [1, 0],    #  x <= 1
              [0, -1],   # -y <= 0  (y >= 0)
              [0, 1]])   #  y <= 1
b = np.array([0, 1, 0, 1])

# Create polytope
P = pc.Polytope(A, b)

print("Polytope created successfully")
print(f"Vertices: {P.vertices}")
print(f"Volume: {P.volume}")

Polytope Operations

We can perform various operations on polytopes such as intersection, union, and Minkowski sum ?

import numpy as np
from scipy.spatial import ConvexHull

# Define two simple polytopes as sets of vertices
poly1_vertices = np.array([[0, 0], [2, 0], [2, 2], [0, 2]])
poly2_vertices = np.array([[1, 1], [3, 1], [3, 3], [1, 3]])

# Create convex hulls
poly1 = ConvexHull(poly1_vertices)
poly2 = ConvexHull(poly2_vertices)

# Function to check if point is inside polytope
def point_in_polytope(point, polytope):
    A = polytope.equations[:, :-1]
    b = -polytope.equations[:, -1]
    return np.all(A.dot(point) <= b + 1e-10)

# Find intersection points (simplified approach)
intersection_points = []
for vertex in poly1_vertices:
    if point_in_polytope(vertex, poly2):
        intersection_points.append(vertex)
for vertex in poly2_vertices:
    if point_in_polytope(vertex, poly1):
        intersection_points.append(vertex)

print("Intersection points:")
print(np.array(intersection_points))
print(f"Polytope 1 area: {poly1.volume:.4f}")
print(f"Polytope 2 area: {poly2.volume:.4f}")

Intersection points:
[[2 2]
 [1 1]]
Polytope 1 area: 4.0000
Polytope 2 area: 4.0000

Properties of Polytopes

Let's examine various properties of polytopes that can be computed ?

import numpy as np
from scipy.spatial import ConvexHull

# Create a regular hexagon
angles = np.linspace(0, 2*np.pi, 7)[:-1]  # 6 points
radius = 2
vertices = np.column_stack([radius * np.cos(angles), 
                           radius * np.sin(angles)])

hull = ConvexHull(vertices)

# Calculate properties
centroid = np.mean(vertices[hull.vertices], axis=0)
perimeter = np.sum([np.linalg.norm(vertices[hull.vertices[i]] - 
                                  vertices[hull.vertices[(i+1) % len(hull.vertices)]]) 
                    for i in range(len(hull.vertices))])

print("Polytope Properties:")
print(f"Number of vertices: {len(hull.vertices)}")
print(f"Area: {hull.volume:.4f}")
print(f"Perimeter: {perimeter:.4f}")
print(f"Centroid: ({centroid[0]:.4f}, {centroid[1]:.4f})")

Polytope Properties:
Number of vertices: 6
Area: 10.3923
Perimeter: 12.0000
Centroid: (0.0000, 0.0000)

Visualization Techniques

Here's how to create better visualizations of polytopes with custom styling ?

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull

# Create an irregular polytope
np.random.seed(123)
n_points = 8
angles = np.sort(np.random.uniform(0, 2*np.pi, n_points))
radii = np.random.uniform(1, 3, n_points)
vertices = np.column_stack([radii * np.cos(angles), 
                           radii * np.sin(angles)])

hull = ConvexHull(vertices)

# Create visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Plot 1: Basic polytope
ax1.plot(vertices[:, 0], vertices[:, 1], 'ro', markersize=8, label='All points')
hull_vertices = vertices[hull.vertices]
hull_x = np.append(hull_vertices[:, 0], hull_vertices[0, 0])
hull_y = np.append(hull_vertices[:, 1], hull_vertices[0, 1])
ax1.plot(hull_x, hull_y, 'b-', linewidth=2, label='Polytope boundary')
ax1.fill(hull_x, hull_y, 'lightblue', alpha=0.3)
ax1.set_title('Irregular Polytope')
ax1.grid(True)
ax1.legend()
ax1.axis('equal')

# Plot 2: With annotations
ax2.plot(vertices[:, 0], vertices[:, 1], 'ro', markersize=8)
ax2.plot(hull_x, hull_y, 'b-', linewidth=2)
ax2.fill(hull_x, hull_y, 'lightgreen', alpha=0.3)
# Annotate vertices
for i, vertex in enumerate(hull_vertices):
    ax2.annotate(f'V{i}', (vertex[0], vertex[1]), 
                xytext=(5, 5), textcoords='offset points')
ax2.set_title('Annotated Polytope')
ax2.grid(True)
ax2.axis('equal')

plt.tight_layout()
plt.show()

print(f"Polytope created with {len(hull.vertices)} vertices")

Polytope created with 8 vertices

Conclusion

Polytopes are fundamental geometric objects that can be effectively manipulated in Python using libraries like NumPy, SciPy, and specialized polytope packages. The ConvexHull class provides an excellent starting point for basic polytope operations, while advanced libraries offer more sophisticated features for complex geometric computations.