Random walks and electric networks

A walker moves through a hall at random, one step at a time to a neighboring point of a grid, and turns back at the walls. How likely is it to leave by the front door rather than by the back door? Peter Doyle and Laurie Snell describe the answer in their book Random Walks and Electric Networks. Put a resistor of 1 Ω on every step the walker can take, hold the front door at 1 V and ground the back door. The voltage at a point is then the chance that a walker starting there reaches the front door first.

Both follow the same rule. A walker’s chance is the average of the chances at the points it can step to, since it takes each of those steps equally often. The voltage is the average of the voltages at those points, since by Kirchhoff’s current law the currents through the equal resistors around a point add up to zero. Both are 1 at the front door and 0 at the back door, and one function on the grid has these values and follows the rule, so the two are the same.

The example builds the network from a planex2d floor plan. It checks the voltage in a corridor, where the chance is known exactly, and in a hall with a pillar and a partition wall against thousands of simulated walkers.

This page is part of a collection of field and simulation examples using planex2d.

Open this example in the playground, the interactive code editor in the browser, to change the code and run it.

import math

import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve

The network

The grid points with whole-number coordinates inside the floor plan are the points the walker can stand on. Contains with boundary=False keeps those inside and drops those on the walls. Each point has four steps, to its neighbors on the right, left, above and below. A step is kept if the walk between the two points stays inside the floor plan, which is checked at three points along it, so a step across a pillar or a corner is dropped. A step that ends on a door leads out of the hall. The table steps holds for every point the index of each neighbor, −1 for a step into a wall, and −2 for the first door, −3 for the second:

def network(plan, doors):
    """the grid points inside the plan, and for each its four steps: the next point's index, -1 into a wall, or
    -2 - d out through door d"""
    (x0, y0), (x1, y1) = plan.bounding_box
    grid = np.array([(x, y) for x in range(math.ceil(x0), math.floor(x1) + 1)
                     for y in range(math.ceil(y0), math.floor(y1) + 1)], dtype=float)
    points = grid[plan.Contains(grid, boundary=False)]
    index = {tuple(p): k for k, p in enumerate(points)}
    steps = np.full((len(points), 4), -1)
    for k, p in enumerate(points):
        for m, d in enumerate(np.array([(1, 0), (-1, 0), (0, 1), (0, -1)])):
            if not plan.Contains(p + np.outer((0.25, 0.5, 0.75), d), boundary=False).all():
                continue
            q = tuple(p + d)
            if q in index:
                steps[k, m] = index[q]
            else:
                steps[k, m] = next((-2 - n for n, door in enumerate(doors) if door.Distance(q) < 1e-9), -1)
    return points, steps

Kirchhoff’s current law at a point says that its voltage times the number of its steps equals the sum of the voltages at the ends of these steps, where a step through a door ends at that door’s voltage. Over all points, this is one sparse linear system:

def voltage(steps, volts):
    """the voltage at every point, with door d held at volts[d] and a resistor of 1 ohm on every step"""
    n = len(steps)
    rows, cols = np.nonzero(steps >= 0)
    a = sp.csr_matrix((-np.ones(len(rows)), (rows, steps[rows, cols])), shape=(n, n))
    a = a + sp.diags((steps != -1).sum(axis=1).astype(float))
    b = sum(v * (steps == -2 - d).sum(axis=1) for d, v in enumerate(volts))
    return spsolve(a.tocsc(), b.astype(float))

A corridor

A corridor of 40 × 4 with a door across each end has a known answer. Since the doors span its whole width, the chance depends only on how far along the walker is, and a walker at \(x\) leaves by the door at \(x = 0\) with the chance \(1 - x/40\). It is the gambler’s ruin, in which a player with \(x\) coins in a fair game, one coin won or lost per round, goes broke before owning 40 with the same chance.

corridor = px.QuickSketch().Rectangle(40, 4).Face()
points, steps = network(corridor, [px.Line((0, 0), (0, 4)), px.Line((40, 0), (40, 4))])
v = voltage(steps, [1, 0])
print(len(points), "points | largest difference from 1 - x/40", np.abs(v - (1 - points[:, 0] / 40)).max())

117 points | largest difference from 1 - x/40 5.440092820663267e-15

The hall

The hall is 30 × 20, with a round pillar of radius 3 and a partition wall 1 thick that runs from the north wall down to 6 above the floor. The front door is 4 wide in the west wall, and the back door, also 4 wide, sits in the east wall behind the partition:

hall = (px.QuickSketch().Rectangle(30, 20).Face() - px.Circle((12, 12), 3).Face()
        - px.QuickSketch((20, 6)).Rectangle(1, 14).Face())
front, back = px.Line((0, 8), (0, 12)), px.Line((30, 2), (30, 6))
points, steps = network(hall, [front, back])
print(len(points), "points |", (steps == -2).sum(), "steps out by the front door,", (steps == -3).sum(), "by the back")
# px.Plot(hall)

494 points | 5 steps out by the front door, 5 by the back

With the front door at 1 V and the back door grounded, the voltage falls from the front door toward the back door (Figure 1). It drops across the partition, since a walker just behind it has to go around its lower end to reach the front door, while the back door is close:

v = voltage(steps, [1, 0])
for p in ((19, 15), (22, 15)):
    print("at", p, "| voltage", round(v[int(np.argmin(np.hypot(*(points - p).T)))], 4))
# plt.scatter(points[:, 0], points[:, 1], c=v); plt.show()

at (19, 15) | voltage 0.69
at (22, 15) | voltage 0.1599

Figure 1: The voltage in the hall, with the doors marked

The walkers

The simulation lets 4000 walkers start from one point and records the door each of them leaves by. All walkers take their steps together, as arrays. Each picks one of the steps open to it, and those that step out through a door stop. The share that leaves by the front door scatters around the chance by about \(\sqrt{p(1 - p)/4000}\):

def walkers(steps, start, count, rng):
    """the share of count walkers from the point start that leave by the first door"""
    n = len(steps)
    open_steps = np.full((n, 4), -1)
    for k, row in enumerate(steps):
        kept = row[row != -1]
        open_steps[k, :len(kept)] = kept
    choices = (steps != -1).sum(axis=1)
    here, out = np.full(count, start), np.zeros(count, dtype=int)
    walking = np.arange(count)
    while len(walking):
        k = here[walking]
        nxt = open_steps[k, (rng.random(len(walking)) * choices[k]).astype(int)]
        left = nxt < 0
        out[walking[left]] = nxt[left]
        here[walking[~left]] = nxt[~left]
        walking = walking[~left]
    return np.mean(out == -2)


rng = np.random.default_rng(1)
for start in ((5, 10), (15, 5), (25, 15), (24, 3)):
    k = int(np.argmin(np.hypot(*(points - start).T)))
    share, sigma = walkers(steps, k, 4000, rng), math.sqrt(v[k] * (1 - v[k]) / 4000)
    print("from", start, "| voltage", round(v[k], 4), "| walkers", share, "| off by",
          round(abs(share - v[k]) / sigma, 2), "standard deviations")

from (5, 10) | voltage 0.8498 | walkers 0.854 | off by 0.75 standard deviations
from (15, 5) | voltage 0.5993 | walkers 0.6165 | off by 2.22 standard deviations
from (25, 15) | voltage 0.158 | walkers 0.15775 | off by 0.04 standard deviations
from (24, 3) | voltage 0.2263 | walkers 0.2235 | off by 0.42 standard deviations

Along the row 4 above the floor, from the west wall to the east wall, 2000 walkers from every third point follow the voltage (Figure 2):

row = np.flatnonzero(points[:, 1] == 4)
shares = np.array([walkers(steps, k, 2000, rng) for k in row[::3]])
p = v[row[::3]]
print("largest difference", round(np.abs(shares - p).max(), 4), "| in standard deviations",
      round((np.abs(shares - p) / np.sqrt(p * (1 - p) / 2000)).max(), 2))
# plt.plot(points[row, 0], v[row], points[row[::3], 0], shares, "o"); plt.show()

largest difference 0.0156 | in standard deviations 1.46

Figure 2: The walkers against the voltage

See Inside or outside for Contains on many points at once, and Squaring the square for another resistor network that solves a geometry problem.