Figure 1. Six unit disks packed around a unit disk (in grey).

Consider the problem of packing disks of unit radius around a disk of same radius. The maximum packing is reached for 6 disks (see figure on the left) and the packing density is 1 (there is no space between the disks). The problem is more complex in three dimensions: Gregory and Newton asked back in 17th century how many unit spheres could be packed around a unit sphere. A tentative solution is to align the centers of 12 outer spheres with vertices of an icosahedron (see Figure 2). Some space is left between the spheres and one may wonder whether a 13th sphere can be inserted. Newton guessed, from experimenting with wooden spheres, that this was not possible, but the mathematical proof of this came only in 1953 (Schütte und van der Waerden).

Figure 2. Twelve unit outer spheres (in blue) packed around an inner unit sphere (in red). Space is left. Is that enough for a 13th sphere?

A related question is to determine the maximal radius R of N spheres that we would like to pack around sphere of unit radius. It happens that for N=13 this maximal radius R is <1. Its determination and the determination of a/the optimal configuration is a complex optimization problem. In this TD we examine some approaches inspired by statistical mechanics.

Problem

Considering the projection of the outer spheres onto the inner sphere, one is left with the problem of packing disks of radius r on the unit sphere. Denoting by

Figure 3. N=16 disks of radius r on the unit sphere.

the vectors of norm 1 pointing to the positions of the disks center, the non-overlapping conditions amounts to

Figure 4. N=13 non-overlapping spheres of radius R≃0.43 (corresponding to r=0.3).

Explain the relation

between the radius R of the outer spheres and the radius r of the corresponding disks projected on the inner sphere (see Figure 3).

Naive algorithm

In a naive approach, one samples uniformly N points on the surface of the unit sphere. The sample is rejected until the disks do not overlap.

Explain how one can sample uniform points on the surface of the sphere.

Here is an algorithm. The visualisation is based on the Mayavi package.

importrandom,math, numpy
N =13; r =.35; R =1. / (1./r - 1.)def dist(a,b): returnmath.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2)def allowed(positions):
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])printmin(dists)returnmin(dists)>2*r
def unif_sphere():
sigma =1./math.sqrt(3.)
x =[random.gauss(0,sigma)for i inrange(3)]
norm =math.sqrt(sum(xk**2. for xk in x))return[xk/norm for xk in x]test=Truewhiletest:
positions =[unif_sphere()for j inrange(N)]if allowed(positions):
print positions
test=Falsefrom mayavi.mlabimport points3d,show,figure
figure(bgcolor=(0,0,.1),size=(1000,700))
m =1.+R
x,y,z=[m*pos[0]for pos in positions],[m*pos[1]for pos in positions],[m*pos[2]for pos in positions]
rad=[2*R for pos in positions]
x.append(0);y.append(0);z.append(0);rad.append(2)
points3d(x,y,z, rad, scale_factor=1, resolution=100,transparent=True)
show()

Markov Chain Monte Carlo algorithm

An example Monte-Carlo algorithm, with moves preserving detailed balance with respect to the uniform measure on the sphere, consists in trying to change the position of a random disk. Here is an example (the visualisation code, identical to the previous one, is not recalled for conciseness).

importrandom,math, numpy
N =13; r =.25; R =1./(1./r-1.); tmax =1000; sig =.3;def dist(a,b): returnmath.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2)def allowed(positions):
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])printmin(dists)returnmin(dists)>2*r
def allowed_move(positions,k,newpos):
dists =[dist(positions[l], newpos)for l inrange(k)+range(k+1,N)]returnmin(dists)>2*r
def unif_sphere():
sigma =1./math.sqrt(3.)
x =[random.gauss(0,sigma)for i inrange(3)]
norm =math.sqrt(sum(xk**2. for xk in x))return[xk/norm for xk in x]#initialisationtest=Truewhiletest:
positions =[unif_sphere()for j inrange(N)]if allowed(positions):
print positions
test=False#moves
count =0for t inrange(tmax):
k =random.randint(0,N-1)
newpos =[positions[k][0]+random.gauss(0,sig),positions[k][1]+random.gauss(0,sig),\
positions[k][2]+random.gauss(0,sig)]
norm =math.sqrt(sum(xk**2. for xk in newpos))
newpos =[xk/norm for xk in newpos]if allowed_move(positions,k,newpos):
positions = positions[:k]+[newpos]+positions[k+1:]
count +=1print"Rejection rate: ",1.*count/tmax

Simulated annealing

Figure 5. Simulated-annealing run for 13 disks on the unit sphere. The ﬁnal density is η = 0.791393

Figure 6. N=13 non-overlapping spheres of radius R=0.88 obtained by simulated annealing from the program below. The packing is η = 0.758

During a Markov-chain simulation, disks almost never touch. This suggests that we should slightly swell the disks at certain times during the simulation, by a small fraction ɣ of some maximum possible increase that would still keep the conﬁguration legal. Such a swelling is implemented between long phases of standard Monte-Carlo simulation.

How can you swell the disks without generating overlap?
Here is an example algorithm.

importrandom,math, numpy
N =13;r =.05;R =1./(1./r-1.);tmax =800000;sig =.25;gamma =0.05random.seed(137)def dist(a,b): returnmath.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2)def allowed(positions):
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])returnmin(dists)>2*r
def allowed_move(positions,k,newpos):
dists =[dist(positions[l], newpos)for l inrange(k)+range(k+1,N)]returnmin(dists)>2*r
def unif_sphere():
sigma =1./math.sqrt(3.)
x =[random.gauss(0,sigma)for i inrange(3)]
norm =math.sqrt(sum(xk**2. for xk in x))return[xk/norm for xk in x]#initialisationtest=Truewhiletest:
positions =[unif_sphere()for j inrange(N)]if allowed(positions): test=False#moves
count =0for t inrange(tmax):
k =random.randint(0,N-1)
newpos =[positions[k][0]+random.gauss(0,sig),positions[k][1]+random.gauss(0,sig),\
positions[k][2]+random.gauss(0,sig)]
norm =math.sqrt(sum(xk**2. for xk in newpos))
newpos =[xk/norm for xk in newpos]if allowed_move(positions,k,newpos):
positions = positions[:k]+[newpos]+positions[k+1:]
count+=1if t%100==0:
if t%10000==0:print' acceptance rate = ',1.*count/(1+t)
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])
Upsilon =min(dists)/2.
r= r+gamma*(Upsilon-r)
R =1./(1./r-1.)print'density',1.*N/2.*(1.-math.sqrt(1-r*r))print'final r',r
print'final R',R

Controlling the annealing rate

Figure 7. N=13 non-overlapping spheres of radius R ≃ 0.9164678249, now obtained by controlling the annealing rate.

In order to keep a decent sampling of the phase space, a simple solution is to adapt the variance of the Gaussian distribution of the Monte Carlo tentative moves. Below we divide this variance σ by two any time the acceptance rate of Monte Carlo tentatives becomes lower than a fixed threshold (fixed to 0.2).

This allows a better mixing of the spheres during the Monte Carlo phases between swellings. Doing so allows us to reach a configuration with R=0.9164678249 in a few seconds, very close to the best known configuration.

importrandom,math, numpy
N =13; r =.2; R =1./(1./r-1.); tmax =500000; sig =.25; gamma =0.1random.seed(138)def dist(a,b): returnmath.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2)def allowed(positions):
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])returnmin(dists)>2*r
def allowed_move(positions,k,newpos):
dists =[dist(positions[l], newpos)for l inrange(k)+range(k+1,N)]returnmin(dists)>2*r
def unif_sphere():
sigma =1./math.sqrt(3.)
x =[random.gauss(0,sigma)for i inrange(3)]
norm =math.sqrt(sum(xk**2. for xk in x))return[xk/norm for xk in x]#initialisationtest=Truewhiletest:
positions =[unif_sphere()for j inrange(N)]if allowed(positions): test=False#moves
count =0
count2 =0for t inrange(tmax):
k =random.randint(0,N-1)
newpos =[positions[k][0]+random.gauss(0,sig),positions[k][1]+random.gauss(0,sig),\
positions[k][2]+random.gauss(0,sig)]
norm =math.sqrt(sum(xk**2. for xk in newpos))
newpos =[xk/norm for xk in newpos]if allowed_move(positions,k,newpos):
positions = positions[:k]+[newpos]+positions[k+1:]
count+=1
count2+=1if t%100==0:
accept=count2/100.
if accept<.2: sig = sig/2if t%40000==0:
print' local acceptance rate = ', accept,'\tsigma = ', sig,'\tR = ',1./(1./r-1.)
count2 =0
dists =[]for l inrange(N): dists.extend([dist(positions[k],positions[l])for k inrange(l)])
Upsilon =min(dists)/2.
r= r+gamma*(Upsilon-r)
R =1./(1./r-1.)print'density',1.*N/2.*(1.-math.sqrt(1-r*r))print'final r',r
print'final R',R
from mayavi.mlabimport points3d,show,figure
figure(bgcolor=(0,0,.1),size=(1000,700))
m =1.+R
x,y,z =[m*pos[0]for pos in positions],[m*pos[1]for pos in positions],[m*pos[2]for pos in positions]
rad =[2*R for pos in positions]
x.append(0); y.append(0); z.append(0); rad.append(2)
points3d(x,y,z, rad, scale_factor=1, resolution=100, transparent=True)
show()

Class Session 10: Simulated annealing for disks on the sphere## Table of Contents

## Introduction

Consider the problem of packing disks of unit radius around a disk of same radius. The maximum packing is reached for 6 disks (see figure on the left) and the packing density is 1 (there is no space between the disks). The problem is more complex in three dimensions: Gregory and Newton asked back in 17th century how many unit spheres could be packed around a unit sphere. A tentative solution is to align the centers of 12 outer spheres with vertices of an icosahedron (see Figure 2). Some space is left between the spheres and one may wonder whether a 13th sphere can be inserted. Newton guessed, from experimenting with wooden spheres, that this was not possible, but the mathematical proof of this came only in 1953 (Schütte und van der Waerden).

A related question is to determine the maximal radius

RofNspheres that we would like to pack around sphere of unit radius. It happens that forN=13 this maximal radiusRis <1. Its determination and the determination of a/the optimal configuration is a complex optimization problem. In this TD we examine some approaches inspired by statistical mechanics.## Problem

Considering the projection of the outer spheres onto the inner sphere, one is left with the problem of packing disks of radius r on the unit sphere. Denoting bythe vectors of norm 1 pointing to the positions of the disks center, the non-overlapping conditions amounts to

Explain the relation

between the radius

Rof the outer spheres and the radiusrof the corresponding disks projected on the inner sphere (see Figure 3).## Naive algorithm

In a naive approach, one samples uniformlyNpoints on the surface of the unit sphere. The sample is rejected until the disks do not overlap.## Markov Chain Monte Carlo algorithm

An example Monte-Carlo algorithm, with moves preserving detailed balance with respect to the uniform measure on the sphere, consists in trying to change the position of a random disk. Here is an example (the visualisation code, identical to the previous one, is not recalled for conciseness).## Simulated annealing

During a Markov-chain simulation, disks almost never touch. This suggests that we should slightly swell the disks at certain times during the simulation, by a small fraction

ɣof some maximum possible increase that would still keep the conﬁguration legal. Such a swelling is implemented between long phases of standard Monte-Carlo simulation.How can you swell the disks without generating overlap?

Here is an example algorithm.

## Controlling the annealing rate

In order to keep a decent sampling of the phase space, a simple solution is to adapt the variance of the Gaussian distribution of the Monte Carlo tentative moves. Below we divide this variance

σby two any time the acceptance rate of Monte Carlo tentatives becomes lower than a fixed threshold (fixed to 0.2).This allows a better mixing of the spheres during the Monte Carlo phases between swellings. Doing so allows us to reach a configuration with R=0.9164678249 in a few seconds, very close to the best known configuration.

## References

Das Problem der dreizehn Kugeln, Mathematische Annalen 125, 325 (1953).[Print this page]