Configurations of hard disks in a (periodic) box¶

A question that arises in statistical physics is the structure of randomly packed non-overlapping disks in a region. Suppose you place $N$ non-overlapping disks of radius $\epsilon$ in a box. (For simplicity, let's also assume the box is periodic -- that is disks placed close to an edge "wrap around" to the opposite edge.)

It turns out that if the disks are packed with high enough density, they arrange themselves into a hexagonal structure.

We are going to test this numerically here. To do this, we need to find a way to randomly sample from all possible arrangements of $N$ non-overlapping disks of radius $\epsilon$ uniformly. That is, we need to randomly sample valid configurations in a manner that no configuration is more likely than any other. Once we do this, there is an easy way to measure "how hexagonal" each arrangement is, so we can compute the "average hexagonalness" by performing repeated random trials.

Here we implement a few algorithms that sample valid configurations. The naieve algorithms don't work (we try them below). Metropolis, Rosenbluth, Teller and Teller introduced an algorithm to study exactly this problem, and produce uniformly random samples of valid configurations. This was later generalized by Hastings, and provides a very useful algorithm that that can be used in more generality.

In [1]:
import numpy as np, scipy as sp
from numpy import sqrt, pi, exp, log, floor, ceil, sin, cos, nan, inf
from numpy import linspace, logspace, geomspace, arange
from numpy import empty, zeros, ones, empty_like, zeros_like, ones_like, full, full_like
from numpy import diff, meshgrid, mean, std, argmin, array, eye, amin, amax, fmin, fmax
from numpy.linalg import norm

%matplotlib inline
import matplotlib as mpl
from matplotlib import pylab, mlab, pyplot as plt
from matplotlib.pylab import plot, scatter, contour, xlabel, ylabel, title, legend
from matplotlib.animation import FuncAnimation
plt.rcParams.update({
    'image.cmap': 'coolwarm',
    'animation.html': 'jshtml',
    'animation.embed_limit': 40, # 40 mb
})

import re, multiprocessing as mp, networkx as nx, json, sys, os

from tqdm.notebook import tqdm, trange

rng = np.random.default_rng()
In [2]:
from itertools import combinations
from matplotlib.pylab import Circle, xticks, yticks, subplot


def figsize(sizex, sizey):
    mpl.rcParams['figure.figsize'] = [sizex, sizey]
In [3]:
# Few utility/plotting functions. They use the global variables N, R and ε
# N = number of disks we're using
# ε = radius of the disks
# R = distance below which we should consider disks to be neighbors.


# We're going to describe our configurations as arrays of size (N, 2), so that
# for i ∈ {1, …, N}, x[i] = (x[i,1], x[i, 2]) describes the center of the ith disk
def plot_disks(x, hl=None, ax=None, **kwargs):
    '''Plot a configuration and optinally highlight one disk'''
    if ax is None:
        ax = plt.gca()
    xticks([])
    yticks([])
    Kwargs = dict(ec='k', alpha=.5)
    Kwargs.update(kwargs)
    KwargsHl = Kwargs.copy()
    KwargsHl['fc'] = 'C1'
    for a in x:
        if hl is not None and all(a == hl):
            KWargs = KwargsHl
        else:
            KWargs = Kwargs

        c = a % 1
        ax.add_patch(Circle(c, ε, **KWargs))
        if c[0] < ε:
            ax.add_patch(Circle(c + [1, 0], ε, **KWargs))
        if c[1] < ε:
            ax.add_patch(Circle(c + [0, 1], ε, **KWargs))
        if c[0] > 1-ε:
            ax.add_patch(Circle(c - [1, 0], ε, **KWargs))
        if c[1] > 1-ε:
            ax.add_patch(Circle(c - [0, 1], ε, **KWargs))
        if c[0] < ε and c[1] < ε:
            ax.add_patch(Circle(c + [1, 1], ε, **KWargs))
        if c[0] < ε and c[1] > 1-ε:
            ax.add_patch(Circle(c + [1, -1], ε, **KWargs))
        if c[0] > 1-ε and c[1] > 1-ε:
            ax.add_patch(Circle(c + [-1, -1], ε, **KWargs))
        if c[0] > 1-ε and c[1] < ε:
            ax.add_patch(Circle(c + [-1, 1], ε, **KWargs))


def place_on_grid(type='rectangular'):
    '''Places N disks on a grid (hexagonal or rectangular), and returns their centers'''
    if type == 'rectangular':
        Nx = int(sqrt(N))
    else:
        Nx = int(sqrt(N * 2 / sqrt(3)))
    Ny = int(ceil(N / Nx))
    Vx, dx = linspace(0, 1, num=Nx, endpoint=False, retstep=True)
    Vy = linspace(0, 1, num=Ny, endpoint=False)

    xx, yy = meshgrid(Vx, Vy)
    if type != 'rectangular':
        xx[::2] += dx/2

    return array(list(zip(xx.flat, yy.flat)))[:N]


# The disk density is the total volume of disks. This finds the radius of disks for
# a given density
def disk_radius(η):
    return sqrt(η / (N*pi))


# Compute R -- balls with centers less than R apart are considerd neighbors
# A reasonable guess is that 1/sqrt(N) should be the 'typical distance' between balls.
# Increasing this by a bit (30%) gives reasonable results.
def neighbor_distance(N):
    return 1.3 / sqrt(N)


# Torus distance. Keep in mind that opposite sides of the square [0, 1]^2 are identified.
# So the point (.9, 0) and (.1, 1) are a distance of .2 apart.
# This is a simple way to compute the square of the distance (on the torus) between two points
def dist2(a, b, axis=-1):
    return np.sum((.5 - abs((a - b) % 1 - .5))**2, axis=axis)


def neighbors(x, a, R):
    '''Retuns all ball centers in x which are at most R away from a (and not identical to a)'''
    n = dist2(x, a)
    return x[(n < R**2) * (n > 0)]


def overlapping(x):
    '''Check if the balls are overlapping'''
    for a, b in combinations(x, 2):
        if dist2(a, b) < (2*ε)**2:
            return True
    return False


def overlaps_with(x, i):
    '''Check if any of the balls with centers x[j] overlap with the ball with center x[i]'''
    return any(dist2(x[i != arange(N)], x[i]) < (2*ε)**2)


def separation2(x):
    return amin([dist2(a, b) for a, b in combinations(x, 2)])
In [4]:
N = 30  # For best results, hoose N to be a perfect square, or somethign with nearly equal factors
R = neighbor_distance(N)

ε = disk_radius(.65)

# Start with a uniform grid arrangement. (We're picking a random disk and highlighting
# it and its neighbors)

x0 = place_on_grid('rectangular')
y0 = place_on_grid('hexagonal')

i = rng.integers(N)
figsize(12, 4)
subplot(1, 2, 1)
plot_disks(x0, hl=x0[i])
plot_disks(neighbors(x0, x0[i], R), fc='C2')

subplot(1, 2, 2)
plot_disks(y0, hl=y0[i])
plot_disks(neighbors(y0, y0[i], R), fc='C2')

figsize(6, 4)
# print( F[0], separation2(x0)-(2*ε)**2 )
if overlapping(x0) or overlapping(y0):
    print('Warning: initial configuration overalps')
No description has been provided for this image

Measuring hexagonalness of configurations¶

One function that measures the hexagonal-ness of a configuration is \begin{equation} F(x) = \frac{1}{N} \Bigl| \sum_{j=1}^N \frac{1}{ |\mathcal N_j|} \sum_{k \in \mathcal N_j} e^{6 i \theta_{j,k}} \Bigr| \,, \end{equation} where $\mathcal N_j$ is the set of all neighbors of the $j$-th disk, and $\theta_{j, k}$ is the angle between the line joining the centers of the $j$-th and $k$-th disks, and the horizontal axis. This function is close to $1$ for hexagonal arrangements, and is typically small otherwise.

In [5]:
def torus_diff(a, b):
    return (a - b + .5) % 1 - .5


def hex_measure(x, R):
    '''Returns a measure of how hexagonal the configuration x is, by averaging over 
    neighbors (ie balls whose centers are at most R apart)'''
    F = 0
    for a in x:
        y = neighbors(x, a, R)
        if len(y) > 0:
            F += np.average(exp(6j * np.arctan2(torus_diff(y[:, 1], a[1]), torus_diff(y[:, 0], a[0]))))
    return abs(F) / N
In [6]:
# Do a sanity test. Print hex_measure of x0 and y0; should be small for x0 and close to 1 for y0
print( f'hex_measure(x0)={hex_measure(x0, R)}, hex_measure(y0)={hex_measure(y0, R)}' )

hex_measure(x0)=1.8369701987210294e-16, hex_measure(y0)=0.9966076395956307

Rejection sampling¶

One algorithm to sample valid configurations uniformly is to place $N$ disks uniformly at random. If they overlap, then reject it and try again.

When the disk density $\eta = N\pi \epsilon^2$ is tiny (about $0.07$), the acceptance rate is 1%, and this algorithm is practical. As $\eta$ approaches the values we care about (0.72), the acceptance rate is so small that this algorithm is computationally intractable.

In [7]:
N = 30
η = .07
ε = sqrt(η / (N*pi))

samples = []
n_trials = 10**3
for _ in trange(n_trials):
    x = rng.uniform(0, 1, size=(N, 2))
    if not overlapping(x):
        samples.append(x)

i = rng.integers(N)
figsize(12, 4)
subplot(1, 2, 1)
plot_disks(samples[0])

subplot(1, 2, 2)
plot_disks(samples[1])

figsize( 6, 4)
print(f'Accepted {len(samples)} in {n_trials} (rate {len(samples) / n_trials * 100}%)')

  0%|          | 0/1000 [00:00<?, ?it/s]
Accepted 11 in 1000 (rate 1.0999999999999999%)
No description has been provided for this image

Rejection sampling, by adding one disk at a time¶

  1. Place one disk by choosing its center uniformly at random.
  2. Given $k$ non-overlapping disks, place one more disk by choosing its center uniformly at random.
    1. If the disk overlaps with any existing disk, reject this configuration and try again.
    2. If the disk does not overlap with any existing disk, we now have a configuration of  $k+1$ non-overlapping disks.
  3. Repeat until we get $N$ non-overlapping disks of radius $\epsilon$.

While this seems like an intuitive approach, the algorithm doesn't work.

  1. While this algorithm gives you $N$ randomly packed non-overlapping disks of radius $\epsilon$, it does not guarantee the configurations are sampled uniformly!
  2. The algorithm may also get stuck. For instance, when the disk density $N \pi \epsilon^2 $ is large, for some $k < N$ you may have chosen a configuration for which it is impossible to add one more non-overlapping disk. The above algorithm will now get stuck, and will need to backtrack and throw away one of the existing disks and try again. This becomes extremely expensive, and the cost of producing even one sample using this algorithm is not computationally tractable.
In [8]:
# Rejection by adding one disk at a time
N = 100
η = .6
ε = sqrt(η / (N*pi))

n_tries = 100
x = rng.uniform(0, 1, size=(1, 2))
for n in trange(N-1):
    for _ in range(n_tries):
        p = rng.uniform(0, 1, size=(1, 2))
        y = np.concat((x, p))
        if not overlapping(y):
            x = y
            break
    else:
        print("Couldn't find a non-overlapping disk")
        break

plot_disks(x)
_ = title(f'{len(x)} disks')

  0%|          | 0/99 [00:00<?, ?it/s]
Couldn't find a non-overlapping disk
No description has been provided for this image

Metropolis Hastings Algorithm¶

The Metropolis Hastings algorithm in this context is as follows:

  1. Start with a valid configuration (for instance where the disks are on a rectangular grid).
  2. Move one of the disks slightly; if it gives you a valid configuration, accept it. If not, reject it and try again.

As you can see the above describes a Markov chain on valid configurations. One can prove that the stationary distribution of this chain is the desired distribution we wish to sample from. Now our sampling problem reduces to iterating this chain "many times" till it gets close to the sationary distribution.

The nice thing about this approach is that it applies to many different situations. The downside is that sometimes the convergence to the stationary distribution may happen very vry slowly, so you need to perform a lot of iterations, which is computationally expensive.

In [9]:
def sample(x0, n_iters, F, show_progress=True):
    x = x0.copy()
    if show_progress:
        accepted = 0
        pbar = tqdm(range(n_iters))
    else:
        pbar = range(n_iters)
    for _ in pbar:
        i = rng.integers(N)
        y = x.copy()
        # Sample uniformly from a ball of radius ρ
        ρ = 1/sqrt(N) - 2*ε
        r = sqrt(rng.uniform(0, ρ**2))
        θ = rng.uniform( 0, 2*pi)
        y[i] += [r*cos(θ), r*sin(θ)]
        y[i] %= 1
        if not overlaps_with(y, i):
            x = y
            F.append(hex_measure(x, R))
            if show_progress:
                accepted += 1
                pbar.set_description(f'Accepted {accepted}, F={F[-1]:.2f}', refresh=False)
        else:
            F.append(F[-1])

    return x, F

When η is too large, the disks essentially become locked in place. The Markov chain is not irreducible, and you're not guaranteed convergence to the distribution you're interested in. Also, since we start with a regular arrangment of the disks, the first few iterations of sample won't be "properly random". So when measuring anything, you should throw away the first few iterations (this is called the burn in period).

Parallelizing: There are a few ways to measure an average of $F$ from this algorithm:

  1. Run it for n_iters iterations; throw away the first few, and then average F. This is hard to parallelize.
  2. Run M independent trials for n_iters/M iterations; throw away the first few iterations (from all trials), then average $F$ over all trials and all the remaining iterations. Same amount of work, but you can parallelize the trials!
In [10]:
# Let's start with a rectangular grid, high density and see
# if the typical configuration becomes hexagonal.
# It's usually easier if you choose N to be a perfect square. You also have a longer
# burn in period. I'm plotting all iterations, but only averaging F after the burn in period
N = 100
n_iters = 10**6
R = neighbor_distance(N)
x0 = place_on_grid('rectangular')
ε = disk_radius(.72)
x, F = sample(x0, n_iters, [hex_measure(x0, R)])

figsize(12, 4)
subplot(1, 2, 1)
plot_disks(x)

subplot(1, 2, 2)
plot(F)
title(f'Average $F={np.average(F[20000:])}$')

figsize(6, 4)

  0%|          | 0/1000000 [00:00<?, ?it/s]
No description has been provided for this image
In [11]:
plot_disks(x)
#savefig( 'figures/hard-disks-eta-72.svg', bbox_inches='tight' )
No description has been provided for this image
In [12]:
# Same as above with lower disk density
ε = disk_radius(.5)
x, F = sample(x0, n_iters, [hex_measure(x0, R)])

figsize(12, 4)
subplot(1, 2, 1)
plot_disks(x)

subplot(1, 2, 2)
plot(F)
title(f'Average $F={np.average(F[20000:])}$')

figsize(6, 4)

  0%|          | 0/1000000 [00:00<?, ?it/s]
No description has been provided for this image
In [13]:
plot_disks(x)
#savefig( 'figures/hard-disks-eta-50.svg', bbox_inches='tight' )
No description has been provided for this image

Some animations.¶

These are just to help visualize the sampling.

In [14]:
import matplotlib.animation as animation
plt.rcParams["animation.html"] = "jshtml"
plt.rcParams['animation.embed_limit'] = 40  # 40 mb
plt.rcParams['savefig.bbox'] = 'tight'
plt.rcParams['savefig.pad_inches'] = 0
In [15]:
N = 36
R = neighbor_distance(N)
ε = disk_radius(.5)
x0 = place_on_grid('rectangular')

fig, ax = plt.subplots()

iters_per_frame = N
(x, F) = (x0, [hex_measure(x0, R)])


def animate(i):
    global x, F

    ax.clear()
    if i > 0:
        (x, F) = sample(x, iters_per_frame, F, show_progress=False)
    plot_disks(x, ax=ax)
    ax.set_title(f'average(F) = {np.average(F):.2f} after {i*iters_per_frame} iterations')


ani = animation.FuncAnimation(fig, animate, interval=40, blit=False, frames=trange(200))
plt.close()

ani

  0%|          | 0/200 [00:00<?, ?it/s]
Out[15]:
No description has been provided for this image
No description has been provided for this image
In [16]:
# Higher density. Use more iterations per frame as the disks move less.
ε = disk_radius(.72)

fig, ax = plt.subplots()

iters_per_frame = 15*N
(x, F) = (x0, [hex_measure(x0, R)])


def animate(i):
    global x, F

    ax.clear()
    if i > 0:
        (x, F) = sample(x, iters_per_frame, F, show_progress=False)
    plot_disks(x, ax=ax)
    ax.set_title(f'average(F) = {np.average(F):.2f} after {i*iters_per_frame} iterations')


ani = animation.FuncAnimation(fig, animate, interval=40, blit=False, frames=trange(200))
plt.close()

ani

  0%|          | 0/200 [00:00<?, ?it/s]
Out[16]:
No description has been provided for this image
No description has been provided for this image
In [ ]: