Graph Algorithms¶
Sparse matrices for graph representation and analysis.
Adjacency Matrix¶
import numpy as np
from scipy import sparse
def main():
# Graph edges: (from, to)
edges = [(0, 1), (0, 2), (1, 2), (2, 3)]
n_nodes = 4
row = [e[0] for e in edges]
col = [e[1] for e in edges]
data = np.ones(len(edges))
# Undirected: symmetrize
A = sparse.csr_matrix((data, (row, col)), shape=(n_nodes, n_nodes))
A = A + A.T
print("Adjacency matrix:")
print(A.toarray())
if __name__ == "__main__":
main()
Graph Laplacian¶
import numpy as np
from scipy import sparse
from scipy.sparse import linalg as splinalg
def main():
A = sparse.csr_matrix([[0, 1, 1, 0],
[1, 0, 1, 0],
[1, 1, 0, 1],
[0, 0, 1, 0]])
# Degree matrix
D = sparse.diags(np.array(A.sum(axis=1)).flatten())
# Laplacian
L = D - A
# Smallest eigenvalues (connectivity)
vals, vecs = splinalg.eigsh(L.astype(float), k=2, which='SM')
print(f"Smallest eigenvalues: {vals}")
if __name__ == "__main__":
main()
Graph connected iff second smallest Laplacian eigenvalue > 0.
Runnable Example: delaunay_triangulation.py¶
"""
Delaunay Triangulation and Spatial Data Structures - A Practical Guide
Exploring scipy.spatial.Delaunay, ConvexHull, Voronoi concepts, and their
applications in computational geometry. Run this file to learn spatial structures!
"""
import numpy as np
from scipy import spatial
if __name__ == "__main__":
print("=" * 70)
print("DELAUNAY TRIANGULATION AND SPATIAL DATA STRUCTURES")
print("=" * 70)
# ============================================================================
# EXAMPLE 1: Basic Delaunay Triangulation in 2D
# ============================================================================
print("\n1. BASIC DELAUNAY TRIANGULATION IN 2D")
print("-" * 70)
# Create a set of 2D points
np.random.seed(42)
points_2d = np.array([
[0, 0],
[1, 0],
[1, 1],
[0, 1],
[0.5, 0.5],
])
print("Points for triangulation:")
for i, point in enumerate(points_2d):
print(f" Point {i}: {point}")
# Compute Delaunay triangulation
delaunay_2d = spatial.Delaunay(points_2d)
print("\nDelaunay triangulation computed!")
print(f"Number of simplices (triangles): {len(delaunay_2d.simplices)}")
print(f"Each simplex is a triangle with 3 vertex indices\n")
# Display triangles
print("Triangles (simplices):")
for i, triangle in enumerate(delaunay_2d.simplices):
v0, v1, v2 = triangle
p0, p1, p2 = points_2d[v0], points_2d[v1], points_2d[v2]
print(f" Triangle {i}: vertices {triangle}")
print(f" Coordinates: {p0}, {p1}, {p2}")
# ============================================================================
# EXAMPLE 2: Accessing Delaunay Data Structures
# ============================================================================
print("\n2. ACCESSING DELAUNAY STRUCTURE INFORMATION")
print("-" * 70)
# simplices: array of vertex indices for each triangle
simplices = delaunay_2d.simplices
print(f"Simplices shape: {simplices.shape}")
print(f"(rows=triangles, cols=vertices per triangle)")
print(f"Simplices:\n{simplices}")
# points: the original input points
print(f"\nInput points shape: {delaunay_2d.points.shape}")
print(f"Points stored in triangulation:\n{delaunay_2d.points}")
# neighbors: for each simplex, which neighboring simplices share edges
neighbors = delaunay_2d.neighbors
print(f"\nNeighbors array shape: {neighbors.shape}")
print(f"neighbors[i, j] = simplex index sharing edge j of triangle i")
print(f"(value -1 means no neighbor, edge is on convex hull)")
print(f"Neighbors:\n{neighbors}")
# vertices: same as simplices (alternative name)
print(f"\nVertices shape: {delaunay_2d.vertices.shape}")
# ============================================================================
# EXAMPLE 3: Understanding the Empty Circumcircle Property
# ============================================================================
print("\n3. DELAUNAY PROPERTY: EMPTY CIRCUMCIRCLE")
print("-" * 70)
print("The Delaunay triangulation has a key property:")
print("The circumcircle of each triangle contains no other points.\n")
def circumcircle(p0, p1, p2):
"""
Compute center and radius of circumcircle for triangle.
Parameters
----------
p0, p1, p2 : array-like
Vertices of the triangle (2D points)
Returns
-------
center : ndarray
Center of circumcircle
radius : float
Radius of circumcircle
"""
# Convert to homogeneous coordinates for calculation
ax, ay = p0
bx, by = p1
cx, cy = p2
d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
if abs(d) < 1e-10:
# Degenerate case
return np.array([0, 0]), 0
ux = ((ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)) / d
uy = ((ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)) / d
center = np.array([ux, uy])
radius = np.linalg.norm(p0 - center)
return center, radius
# Check circumcircle property for first triangle
triangle_0 = delaunay_2d.simplices[0]
v0, v1, v2 = triangle_0
p0, p1, p2 = points_2d[v0], points_2d[v1], points_2d[v2]
center, radius = circumcircle(p0, p1, p2)
print(f"First triangle vertices: {triangle_0}")
print(f" Coordinates: {p0}, {p1}, {p2}")
print(f" Circumcircle center: {center}")
print(f" Circumcircle radius: {radius:.6f}\n")
# Check distances of all points to this circumcircle
print("Distances from other points to circumcircle:")
for i, point in enumerate(points_2d):
dist_to_center = np.linalg.norm(point - center)
if i not in triangle_0:
status = "OUTSIDE" if dist_to_center > radius + 1e-6 else "INSIDE/ON"
print(f" Point {i}: distance to center = {dist_to_center:.6f}, {status}")
# ============================================================================
# EXAMPLE 4: Convex Hull - Boundary of Point Set
# ============================================================================
print("\n4. CONVEX HULL: BOUNDARY OF POINT SET")
print("-" * 70)
print("The convex hull is the smallest convex polygon containing all points.\n")
# Compute convex hull
hull = spatial.ConvexHull(points_2d)
print(f"Convex hull computed!")
print(f"Number of vertices on hull: {len(hull.vertices)}")
print(f"Vertex indices: {hull.vertices}")
# Hull points
print(f"\nHull vertices (coordinates):")
for idx in hull.vertices:
print(f" Point {idx}: {points_2d[idx]}")
# Hull simplices (in 2D, these are edges)
print(f"\nHull simplices (edges in 2D):")
for i, edge in enumerate(hull.simplices):
v0, v1 = edge
p0, p1 = points_2d[v0], points_2d[v1]
print(f" Edge {i}: points {v0} -- {v1}")
# Area and volume
print(f"\nHull properties:")
print(f" Area: {hull.volume:.6f}") # In 2D, 'volume' is actually area
# ============================================================================
# EXAMPLE 5: Finding Hull Edges from Delaunay Triangulation
# ============================================================================
print("\n5. IDENTIFYING BOUNDARY EDGES IN DELAUNAY")
print("-" * 70)
print("In Delaunay, boundary edges have neighbors = -1\n")
boundary_edges = []
for i, neighbor_row in enumerate(delaunay_2d.neighbors):
for edge_idx, neighbor_simplex in enumerate(neighbor_row):
if neighbor_simplex == -1: # No neighbor = boundary edge
# Find vertices of this edge
# Each edge is opposite to a vertex
edge_vertices = np.delete(delaunay_2d.simplices[i], edge_idx)
boundary_edges.append(edge_vertices)
print(f"Found {len(boundary_edges)} boundary edges:")
for i, edge in enumerate(boundary_edges):
v0, v1 = edge
print(f" Edge {i}: points {v0} -- {v1}, coordinates {points_2d[v0]} -- {points_2d[v1]}")
# ============================================================================
# EXAMPLE 6: Point Location in Triangulation
# ============================================================================
print("\n6. POINT LOCATION: FINDING WHICH TRIANGLE CONTAINS A POINT")
print("-" * 70)
print("Query: Which triangle contains a given point?\n")
test_points = np.array([
[0.5, 0.5], # Inside the triangulation
[2.0, 2.0], # Outside
[0.1, 0.1], # Inside
])
print("Test points and their containing triangles:")
for test_point in test_points:
simplex_idx = delaunay_2d.find_simplex(test_point)
if simplex_idx >= 0:
triangle = delaunay_2d.simplices[simplex_idx]
print(f" Point {test_point}: in triangle {simplex_idx}, vertices {triangle}")
else:
print(f" Point {test_point}: OUTSIDE triangulation")
# ============================================================================
# EXAMPLE 7: 3D Delaunay Triangulation with Random Data
# ============================================================================
print("\n7. DELAUNAY TRIANGULATION IN 3D")
print("-" * 70)
print("Extend to 3D: triangles become tetrahedra!\n")
# Generate random 3D points
np.random.seed(42)
points_3d = np.random.rand(20, 3) # 20 random points in unit cube
# Compute 3D Delaunay
delaunay_3d = spatial.Delaunay(points_3d)
print(f"3D Delaunay triangulation:")
print(f" Number of points: {delaunay_3d.points.shape[0]}")
print(f" Number of tetrahedra: {len(delaunay_3d.simplices)}")
print(f" Tetrahedron size: {delaunay_3d.simplices.shape[1]} vertices")
# Show first few tetrahedra
print(f"\nFirst 3 tetrahedra:")
for i in range(min(3, len(delaunay_3d.simplices))):
tet = delaunay_3d.simplices[i]
print(f" Tetrahedron {i}: vertices {tet}")
# ============================================================================
# EXAMPLE 8: Volume Calculation in 3D Delaunay
# ============================================================================
print("\n8. COMPUTING TETRAHEDRON VOLUMES")
print("-" * 70)
def tetrahedron_volume(p0, p1, p2, p3):
"""
Compute volume of tetrahedron with vertices p0, p1, p2, p3.
Volume = |det([p1-p0, p2-p0, p3-p0])| / 6
"""
matrix = np.array([
p1 - p0,
p2 - p0,
p3 - p0
]).T
volume = abs(np.linalg.det(matrix)) / 6.0
return volume
# Compute volumes of first 5 tetrahedra
print("Volumes of first 5 tetrahedra:\n")
total_volume = 0
for i in range(min(5, len(delaunay_3d.simplices))):
tet_vertices = delaunay_3d.simplices[i]
pts = delaunay_3d.points[tet_vertices]
vol = tetrahedron_volume(pts[0], pts[1], pts[2], pts[3])
total_volume += vol
print(f" Tetrahedron {i}: volume = {vol:.8f}")
print(f"\nTotal volume of first 5 tetrahedra: {total_volume:.8f}")
# ============================================================================
# EXAMPLE 9: Voronoi Diagram Concept
# ============================================================================
print("\n9. VORONOI DIAGRAM: DUAL OF DELAUNAY")
print("-" * 70)
print("Voronoi diagram is dual to Delaunay:")
print("- Each Voronoi cell is the region closer to one point than any other")
print("- Voronoi vertices are circumcenters of Delaunay triangles")
print("- Voronoi edges connect these circumcenters\n")
# Compute Voronoi diagram
vor = spatial.Voronoi(points_2d)
print(f"Voronoi diagram for 2D points:")
print(f" Number of Voronoi vertices: {vor.vertices.shape[0]}")
print(f" Number of Voronoi regions: {len(vor.regions)}")
print(f"\nVoronoi vertices (first 3):")
for i in range(min(3, vor.vertices.shape[0])):
print(f" Vertex {i}: {vor.vertices[i]}")
print(f"\nFirst 3 Voronoi regions:")
for i in range(min(3, len(vor.regions))):
region = vor.regions[i]
print(f" Point {i} region: vertices {region}")
# ============================================================================
# EXAMPLE 10: Practical Application - Mesh Generation
# ============================================================================
print("\n10. PRACTICAL APPLICATION: MESH GENERATION")
print("-" * 70)
print("Delaunay triangulation is used for mesh generation in FEA/FEM.\n")
# Create a grid of points
x = np.linspace(-1, 1, 5)
y = np.linspace(-1, 1, 5)
xx, yy = np.meshgrid(x, y)
grid_points = np.column_stack([xx.ravel(), yy.ravel()])
# Create Delaunay mesh
mesh = spatial.Delaunay(grid_points)
print(f"Mesh from 5x5 grid:")
print(f" Points: {mesh.points.shape[0]}")
print(f" Triangles: {len(mesh.simplices)}")
# Show triangle quality metrics
def triangle_aspect_ratio(p0, p1, p2):
"""Compute aspect ratio of triangle (1.0 = equilateral)"""
side_lengths = [
np.linalg.norm(p1 - p0),
np.linalg.norm(p2 - p1),
np.linalg.norm(p0 - p2)
]
return max(side_lengths) / min(side_lengths)
aspect_ratios = []
for triangle in mesh.simplices:
p0, p1, p2 = mesh.points[triangle]
ar = triangle_aspect_ratio(p0, p1, p2)
aspect_ratios.append(ar)
aspect_ratios = np.array(aspect_ratios)
print(f"\nTriangle aspect ratios:")
print(f" Min: {aspect_ratios.min():.4f} (most regular)")
print(f" Max: {aspect_ratios.max():.4f} (most irregular)")
print(f" Mean: {aspect_ratios.mean():.4f}")
print(f" (1.0 = perfect equilateral, higher = worse)")
print("\n" + "=" * 70)
print("Key takeaways:")
print("- Delaunay triangulation: optimal triangulation maximizing angles")
print("- Empty circumcircle property: no points inside triangle's circumcircle")
print("- ConvexHull: boundary of point set, useful for collision detection")
print("- 2D: simplices are triangles, 3D: simplices are tetrahedra")
print("- find_simplex(): locate point within triangulation")
print("- Voronoi: dual structure, useful for spatial analysis")
print("- Applications: mesh generation, interpolation, spatial indexing")
print("=" * 70)