The Mapper Algorithm

Mapper is a topological algorithm that constructs a simplified “skeleton” of high-dimensional data by pulling back a cover of the range of a filter function. Applied to cities, it reveals the branching structure, loops, and central nodes of urban form—a topological X-ray of metropolitan organization.

1. The Mapper Construction

Given a point cloud $X$ and a continuous function $f: X \to \mathbb{R}$ (the filter function), Mapper constructs a simplicial complex that captures the topology of $X$ through four steps:

Step 1: Choose a filter function. The filter maps each data point to a real value. For urban analysis, natural choices include:

Distance from the Central Business District (CBD).
Population density or land value.
Eccentricity: the maximum geodesic distance to any other point.
A principal component of socioeconomic data.

Step 2: Cover the range. Choose an open cover $\mathcal{U} = \{U_1, U_2, \dots, U_n\}$ of $f(X) \subseteq \mathbb{R}$. Typically, overlapping intervals of width $w$ with overlap fraction $p$:

$$U_i = [a_i, a_i + w], \quad a_{i+1} = a_i + w(1 - p)$$

Step 3: Pull back and cluster. For each cover element $U_i$, consider the preimage $f^{-1}(U_i) \cap X$—all data points whose filter value falls in $U_i$. Apply a clustering algorithm (e.g., single-linkage, DBSCAN) to partition each preimage into clusters:

$$f^{-1}(U_i) \cap X = C_{i,1} \sqcup C_{i,2} \sqcup \cdots \sqcup C_{i,k_i}$$

Step 4: Build the nerve. Create a vertex for each cluster $C_{i,j}$. Connect two vertices with an edge if their clusters share data points (which can happen when cover elements overlap). More generally, form a $k$-simplex whenever $k+1$ clusters have a common point.

The resulting simplicial complex is the Mapper graph—a topological skeleton that captures the essential shape of the data as seen through the filter function.

2. Mathematical Foundation: The Nerve Theorem

Mapper is grounded in the Nerve Theorem, which guarantees that the nerve of a “good” cover captures the topology of the underlying space:

$$\text{If } \mathcal{U} = \{U_i\} \text{ is a good cover of } X, \text{ then } \mathcal{N}(\mathcal{U}) \simeq X$$

where $\mathcal{N}(\mathcal{U})$ is the nerve and $\simeq$ denotes homotopy equivalence. A cover is “good” if all finite intersections $U_{i_1} \cap \cdots \cap U_{i_k}$ are either empty or contractible.

In Mapper, we work with the pullback cover $\{f^{-1}(U_i)\}$. The clustering step refines each preimage into connected components, ensuring the cover elements are as close to contractible as the data allows. The quality of the Mapper output depends on:

Filter function: Should capture meaningful variation in the data. Distance from CBD captures the monocentric structure.
Cover resolution: Number and overlap of intervals. More intervals give finer resolution but more complex output.
Clustering parameters: Control how aggressively points are grouped. Affects the granularity of the skeleton.

3. Urban Interpretation of Mapper Features

The Mapper graph of a city reveals its organisational skeleton:

Linear chains: Radial corridors extending from the centre. Each node is a distance band from the CBD, and the linear structure shows how the city extends along major axes.
Branch points: Where the city splits into distinct corridors or sectors. A node with degree > 2 indicates a fork in the urban structure.
Loops: Ring roads, orbital connections, or areas where two corridors reconverge. Loops in the Mapper graph correspond to genuine loops in the urban fabric.
Disconnected components: Satellite towns or isolated developments with no continuous urban connection to the main city at the chosen scale.
Hub nodes (large clusters): Nodes containing many data points represent dense, well-connected areas. The CBD typically appears as the largest hub.

4. Mapper with Distance-from-CBD Filter

We implement the full Mapper algorithm for a synthetic city with a central business district, radial corridors, and a ring road. The filter function is Euclidean distance from the CBD.

Mapper Algorithm with Distance-from-CBD Filter

Python

script.py260 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform
from collections import defaultdict

np.random.seed(42)

# ============================================
# Generate Synthetic City
# ============================================
def generate_city(n=300):
    points = []

# CBD core (dense cluster at origin)
    cbd = np.random.randn(40, 2) * 0.5
    points.append(cbd)

# Radial corridors (4 main axes)
    for angle in [0, np.pi/2, np.pi, 3*np.pi/2]:
        for r in np.linspace(1, 6, 15):
            n_pts = max(2, int(8 - r))
            corridor = np.random.randn(n_pts, 2) * 0.3
            corridor[:, 0] += r * np.cos(angle)
            corridor[:, 1] += r * np.sin(angle)
            points.append(corridor)

# Ring road (circle at r=4)
    n_ring = 50
    angles_ring = np.linspace(0, 2*np.pi, n_ring, endpoint=False)
    ring = np.column_stack([4 * np.cos(angles_ring), 4 * np.sin(angles_ring)])
    ring += np.random.randn(n_ring, 2) * 0.2
    points.append(ring)

# Suburban scatter
    n_suburb = 30
    angles_sub = np.random.rand(n_suburb) * 2 * np.pi
    radii_sub = 5 + np.random.rand(n_suburb) * 3
    suburbs = np.column_stack([radii_sub * np.cos(angles_sub),
                                radii_sub * np.sin(angles_sub)])
    suburbs += np.random.randn(n_suburb, 2) * 0.4
    points.append(suburbs)

return np.vstack(points)

city = generate_city()
n_pts = len(city)
print(f"City point cloud: {n_pts} points")

# ============================================
# Mapper Algorithm
# ============================================

# Step 1: Filter function = distance from CBD (origin)
cbd_location = np.array([0.0, 0.0])
filter_values = np.linalg.norm(city - cbd_location, axis=1)

# Step 2: Cover the range with overlapping intervals
n_intervals = 12
overlap = 0.3
f_min, f_max = filter_values.min(), filter_values.max()
f_range = f_max - f_min
width = f_range / (n_intervals - (n_intervals - 1) * overlap)
step = width * (1 - overlap)

intervals = []
for i in range(n_intervals):
    lo = f_min + i * step
    hi = lo + width
    intervals.append((lo, hi))

print(f"\nFilter range: [{f_min:.2f}, {f_max:.2f}]")
print(f"Intervals: {n_intervals}, width: {width:.2f}, overlap: {overlap}")

# Step 3: Pull back and cluster (single-linkage with distance threshold)
cluster_threshold = 1.5  # distance threshold for clustering
dist_matrix = squareform(pdist(city))

mapper_nodes = []  # list of (interval_idx, cluster_points)
node_id = 0

for i, (lo, hi) in enumerate(intervals):
    # Points in this interval
    mask = (filter_values >= lo) & (filter_values <= hi)
    indices = np.where(mask)[0]

if len(indices) == 0:
        continue

# Single-linkage clustering
    visited = set()
    clusters = []

for idx in indices:
        if idx in visited:
            continue
        # BFS to find cluster
        cluster = []
        queue = [idx]
        while queue:
            node = queue.pop(0)
            if node in visited:
                continue
            visited.add(node)
            cluster.append(node)
            for other in indices:
                if other not in visited and dist_matrix[node, other] < cluster_threshold:
                    queue.append(other)
        if cluster:
            clusters.append(set(cluster))

for cluster_pts in clusters:
        mapper_nodes.append({
            'id': node_id,
            'interval': i,
            'points': cluster_pts,
            'size': len(cluster_pts),
            'center': np.mean(city[list(cluster_pts)], axis=0),
            'mean_filter': np.mean(filter_values[list(cluster_pts)])
        })
        node_id += 1

print(f"\nMapper nodes: {len(mapper_nodes)}")

# Step 4: Build nerve (connect nodes sharing points)
mapper_edges = []
for a in range(len(mapper_nodes)):
    for b in range(a + 1, len(mapper_nodes)):
        shared = mapper_nodes[a]['points'] & mapper_nodes[b]['points']
        if shared:
            mapper_edges.append((a, b, len(shared)))

print(f"Mapper edges: {len(mapper_edges)}")

# Compute Betti numbers of Mapper graph
# Union-Find for beta_0
parent = list(range(len(mapper_nodes)))
rank_uf = [0] * len(mapper_nodes)

def find(x):
    while parent[x] != x:
        parent[x] = parent[parent[x]]
        x = parent[x]
    return x

def union(a, b):
    ra, rb = find(a), find(b)
    if ra == rb: return False
    if rank_uf[ra] < rank_uf[rb]: ra, rb = rb, ra
    parent[rb] = ra
    if rank_uf[ra] == rank_uf[rb]: rank_uf[ra] += 1
    return True

for a, b, _ in mapper_edges:
    union(a, b)

n_components = len(set(find(i) for i in range(len(mapper_nodes))))
n_v = len(mapper_nodes)
n_e = len(mapper_edges)
beta_0 = n_components
beta_1 = n_e - n_v + beta_0

print(f"\n=== Mapper Graph Topology ===")
print(f"Vertices: {n_v}")
print(f"Edges: {n_e}")
print(f"beta_0 (components): {beta_0}")
print(f"beta_1 (loops): {beta_1}")

# Degree distribution
degrees = defaultdict(int)
for a, b, _ in mapper_edges:
    degrees[a] += 1
    degrees[b] += 1

branch_points = sum(1 for d in degrees.values() if d > 2)
endpoints = sum(1 for v in range(n_v) if degrees.get(v, 0) <= 1)
print(f"Branch points (degree > 2): {branch_points}")
print(f"Endpoints (degree <= 1): {endpoints}")

# --- Plot ---
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# Point cloud colored by filter value
ax = axes[0]
scatter = ax.scatter(city[:, 0], city[:, 1], c=filter_values, cmap='viridis',
                     s=10, alpha=0.6)
ax.plot(0, 0, 'r*', markersize=20, label='CBD')
plt.colorbar(scatter, ax=ax, label='Distance from CBD')
ax.set_title('City Point Cloud', color='white', fontweight='bold', fontsize=13)
ax.set_xlabel('x (km)', color='white')
ax.set_ylabel('y (km)', color='white')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_aspect('equal')
ax.grid(True, alpha=0.1)

# Mapper graph in geographic space
ax = axes[1]
for a, b, w in mapper_edges:
    ca = mapper_nodes[a]['center']
    cb = mapper_nodes[b]['center']
    ax.plot([ca[0], cb[0]], [ca[1], cb[1]], '-', color='#10b981', alpha=0.5, linewidth=1.5)

node_sizes = [n['size'] * 3 for n in mapper_nodes]
node_colors = [n['mean_filter'] for n in mapper_nodes]
centers = np.array([n['center'] for n in mapper_nodes])
scatter2 = ax.scatter(centers[:, 0], centers[:, 1], c=node_colors, cmap='viridis',
                      s=node_sizes, alpha=0.8, edgecolors='white', linewidth=0.5, zorder=5)

# Highlight branch points
for v in range(n_v):
    if degrees.get(v, 0) > 2:
        c = mapper_nodes[v]['center']
        ax.scatter(c[0], c[1], s=200, facecolors='none', edgecolors='#ef4444',
                   linewidth=2, zorder=6)

ax.plot(0, 0, 'r*', markersize=20, zorder=7)
ax.set_title('Mapper Graph (geographic)', color='white', fontweight='bold', fontsize=13)
ax.set_xlabel('x (km)', color='white')
ax.set_ylabel('y (km)', color='white')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_aspect('equal')
ax.grid(True, alpha=0.1)

# Mapper graph abstract layout (by filter value)
ax = axes[2]
for a, b, w in mapper_edges:
    fa = mapper_nodes[a]['mean_filter']
    fb = mapper_nodes[b]['mean_filter']
    ya = mapper_nodes[a]['center'][1]
    yb = mapper_nodes[b]['center'][1]
    ax.plot([fa, fb], [ya, yb], '-', color='#10b981', alpha=0.5, linewidth=1.5)

ax.scatter([n['mean_filter'] for n in mapper_nodes],
           [n['center'][1] for n in mapper_nodes],
           c=node_colors, cmap='viridis', s=node_sizes, alpha=0.8,
           edgecolors='white', linewidth=0.5, zorder=5)

ax.set_xlabel('Distance from CBD', color='white', fontsize=12)
ax.set_ylabel('Lateral position (km)', color='white', fontsize=12)
ax.set_title('Mapper Graph (filter space)', color='white', fontweight='bold', fontsize=13)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.1)

fig.patch.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\nPlot saved: Mapper algorithm applied to synthetic city")
print("\nInterpretation:")
print("- Central hub = CBD (large node at distance 0)")
print("- Branches = radial corridors")
print("- Loop at distance ~4 = ring road")
print("- Red circles = branch points where corridors fork")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Union-Find for Betti Numbers

Computing $\beta_0$ (connected components) is efficiently done via the Union-Find (disjoint set) data structure, which supports near-constant-time union and find operations with path compression and union by rank.

The cycle rank formula then gives $\beta_1$ directly:

$$\beta_1 = E - V + \beta_0$$

where $E$ is the number of edges and $V$ the number of vertices. This formula follows from the Euler characteristic: $\chi = V - E = \beta_0 - \beta_1$.

Below is a Fortran-style implementation in Python that demonstrates the Union-Find algorithm applied to a Mapper graph, computing both Betti numbers efficiently.

Union-Find for Betti Numbers (Fortran-style implementation)

Python

script.py222 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

np.random.seed(42)

# ============================================
# Fortran-style Union-Find Implementation
# ============================================

def union_find_betti(n_vertices, edges):
    """
    Compute beta_0 and beta_1 using Union-Find.

Fortran-style: arrays instead of objects, explicit loops,
    path compression via iterative find.
    """
    # Parent and rank arrays (0-indexed)
    parent = np.arange(n_vertices, dtype=np.int32)
    rnk = np.zeros(n_vertices, dtype=np.int32)

# Find with path compression (iterative)
    def find(x):
        root = x
        while parent[root] != root:
            root = parent[root]
        # Path compression
        while parent[x] != root:
            next_x = parent[x]
            parent[x] = root
            x = next_x
        return root

# Union by rank
    def union(a, b):
        ra, rb = find(a), find(b)
        if ra == rb:
            return False  # already in same set
        if rnk[ra] < rnk[rb]:
            ra, rb = rb, ra
        parent[rb] = ra
        if rnk[ra] == rnk[rb]:
            rnk[ra] += 1
        return True

# Process edges
    n_edges = len(edges)
    n_merges = 0
    n_cycles = 0

for e in range(n_edges):
        u, v = edges[e]
        if union(u, v):
            n_merges += 1
        else:
            n_cycles += 1

# Count components
    components = set()
    for i in range(n_vertices):
        components.add(find(i))
    beta_0 = len(components)

# Cycle rank formula
    beta_1 = n_edges - n_vertices + beta_0

return beta_0, beta_1, n_merges, n_cycles

# ============================================
# Test on Multiple Graph Types
# ============================================

def make_grid(side):
    n = side * side
    edges = []
    for i in range(side):
        for j in range(side):
            v = i * side + j
            if j < side - 1: edges.append((v, v + 1))
            if i < side - 1: edges.append((v, v + side))
    return n, edges

def make_tree(depth, branching=2):
    n = (branching**(depth+1) - 1) // (branching - 1)
    edges = []
    for i in range(n):
        for b in range(branching):
            child = branching * i + b + 1
            if child < n:
                edges.append((i, child))
    return n, edges

def make_cycle_graph(n):
    edges = [(i, (i+1) % n) for i in range(n)]
    return n, edges

def make_complete(n):
    edges = [(i, j) for i in range(n) for j in range(i+1, n)]
    return n, edges

def make_barbell(n_each):
    n = 2 * n_each + 1
    edges = []
    # Left clique
    for i in range(n_each):
        for j in range(i+1, n_each):
            edges.append((i, j))
    # Right clique
    for i in range(n_each + 1, n):
        for j in range(i+1, n):
            edges.append((i, j))
    # Bridge
    edges.append((0, n_each))
    edges.append((n_each, n_each + 1))
    return n, edges

test_graphs = {
    'Grid 5x5': make_grid(5),
    'Grid 10x10': make_grid(10),
    'Tree (d=4, b=3)': make_tree(4, 3),
    'Cycle (n=20)': make_cycle_graph(20),
    'Complete K_8': make_complete(8),
    'Barbell (5+5)': make_barbell(5),
}

print("=" * 65)
print(f"{'Graph':<20} {'V':>5} {'E':>5} {'beta_0':>7} {'beta_1':>7} {'chi':>5}")
print("=" * 65)

results = {}
for name, (n_v, edges) in test_graphs.items():
    b0, b1, merges, cycles = union_find_betti(n_v, edges)
    chi = b0 - b1
    results[name] = (n_v, len(edges), b0, b1, chi)
    print(f"{name:<20} {n_v:>5} {len(edges):>5} {b0:>7} {b1:>7} {chi:>5}")

print("=" * 65)

# Verify: for a tree, beta_1 = 0
print("\n--- Verification ---")
print("Tree beta_1 should be 0:", results['Tree (d=4, b=3)'][3] == 0)
print("Cycle beta_1 should be 1:", results['Cycle (n=20)'][3] == 1)
print("Complete K_8 beta_1 = E - V + 1 =", results['Complete K_8'][1] - results['Complete K_8'][0] + 1)

# ============================================
# Scaling Test: Union-Find Performance
# ============================================
import time

sizes = [100, 500, 1000, 2000, 5000]
times_uf = []

for n_side in [10, 22, 32, 45, 71]:  # approximate sqrt of sizes
    n_v, edges = make_grid(n_side)
    start = time.time()
    for _ in range(10):  # repeat for timing
        b0, b1, _, _ = union_find_betti(n_v, edges)
    elapsed = (time.time() - start) / 10
    times_uf.append((n_v, len(edges), elapsed, b0, b1))
    print(f"Grid {n_side}x{n_side}: V={n_v}, E={len(edges)}, "
          f"beta_0={b0}, beta_1={b1}, time={elapsed*1000:.2f}ms")

# --- Plot ---
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Betti number comparison
ax = axes[0]
names = list(results.keys())
b0_vals = [results[n][2] for n in names]
b1_vals = [results[n][3] for n in names]
x = np.arange(len(names))
w = 0.35
ax.bar(x - w/2, b0_vals, w, color='#10b981', alpha=0.8, label='beta_0')
ax.bar(x + w/2, b1_vals, w, color='#06b6d4', alpha=0.8, label='beta_1')
ax.set_xticks(x)
ax.set_xticklabels([n.replace(' ', '\n') for n in names], fontsize=8, color='white', rotation=30)
ax.set_ylabel('Betti Number', fontsize=12, color='white')
ax.set_title('Betti Numbers by Graph Type', fontsize=13, color='white', fontweight='bold')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2, axis='y')

# Union-Find scaling
ax = axes[1]
v_sizes = [t[0] for t in times_uf]
e_sizes = [t[1] for t in times_uf]
run_times = [t[2] * 1000 for t in times_uf]  # ms
ax.loglog(e_sizes, run_times, 'o-', color='#10b981', markersize=8, linewidth=2)
# Fit line
log_e = np.log(e_sizes)
log_t = np.log(run_times)
slope = np.polyfit(log_e, log_t, 1)[0]
ax.set_xlabel('Number of Edges', fontsize=12, color='white')
ax.set_ylabel('Time (ms)', fontsize=12, color='white')
ax.set_title(f'Union-Find Scaling (slope={slope:.2f})', fontsize=13, color='white', fontweight='bold')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Euler characteristic
ax = axes[2]
chi_vals = [results[n][4] for n in names]
colors_chi = ['#10b981' if c >= 0 else '#ef4444' for c in chi_vals]
ax.barh(range(len(names)), chi_vals, color=colors_chi, alpha=0.8)
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names, fontsize=9, color='white')
ax.set_xlabel('Euler Characteristic (chi = beta_0 - beta_1)', fontsize=11, color='white')
ax.set_title('Euler Characteristic', fontsize=13, color='white', fontweight='bold')
ax.axvline(x=0, color='white', linewidth=0.5, alpha=0.3)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2, axis='x')

fig.patch.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print("\nPlot saved: Union-Find Betti numbers and scaling")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Mapper Parameters and Sensitivity

The Mapper output depends on three main parameters. Understanding their effect is essential for reliable urban analysis:

Number of intervals $n$: Controls resolution. Too few intervals produce a coarse skeleton that misses branching structure. Too many create a complex graph dominated by noise. Typical range: 8-30 for urban datasets.
Overlap fraction $p$: Controls connectivity. Higher overlap means more edges in the Mapper graph (more chances for shared points). Typical range: 20-50%.
Clustering threshold $\varepsilon_c$: Controls granularity within each interval. Smaller threshold means more, smaller clusters (more nodes). Must be calibrated to the spatial scale of the data.

The statistical Mapper framework provides stability results: for sufficiently large datasets, the Mapper graph converges to the Reeb graph of the filter function, which is a topological invariant of the data:

$$\text{Reeb}_f(X) = X / \sim \quad \text{where } x \sim y \iff f(x) = f(y) \text{ and } x, y \text{ connected in } f^{-1}(f(x))$$

The Reeb graph collapses each connected component of a level set to a single point, producing the simplest possible representation of the topology through the filter lens. Mapper approximates this Reeb graph for finite data, with the approximation quality improving as the number of data points grows.

7. The Topological Toolkit for Urban Analysis

Modules 15 has assembled a complete topological toolkit for urban analysis:

Simplicial homology: Computes Betti numbers at a single scale. Best for comparing canonical urban forms.
Persistent homology: Tracks topological features across scales. Reveals multi-scale structure and separates signal from noise.
Mapper: Constructs a topological skeleton through a chosen filter. Reveals branching, loops, and hubs in the urban fabric.

Together, these tools provide a coordinate-free, scale-invariant, and deformation-robust language for describing urban form. Unlike metric measures (street density, block area), topological invariants capture the qualitative structure of a city: whether it has loops, branches, disconnections, or enclosed regions.

Key Takeaway

The Mapper algorithm produces a topological skeleton of the city by pulling back a cover of a filter function's range. With distance-from-CBD as the filter, it reveals radial corridors (linear chains), ring roads (loops), branching points (high-degree nodes), and the hierarchical organisation of the metropolitan area. Combined with Union-Find for efficient Betti number computation, it provides a practical and theoretically grounded approach to understanding urban form.

← Persistent Homology Part 16: Information Theory →