Simulated+annealing+for+disks+on+the+sphere

toc =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).
 * Class Session 10: Simulated annealing for disks on the 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 math \vec x_k \:, \qquad 1\leq k\leq N math the vectors of norm 1 pointing to the positions of the disks center, the non-overlapping conditions amounts to math |\vec x_k-\vec x_l|>2r~,\quad\quad k\neq l math Explain the relation math \displaystyle R=\frac{1}{\frac 1 r-1} math 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. > code format="python" import random, math, numpy N = 13; r = .35; R = 1. / (1./r - 1.) def dist(a,b): return math.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2) def allowed(positions): dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) print min(dists) return min(dists) > 2*r def unif_sphere: sigma = 1./math.sqrt(3.) x = [random.gauss(0,sigma) for i in range(3)] norm = math.sqrt(sum(xk**2. for xk in x)) return [xk/norm for xk in x] test=True while test: positions = [unif_sphere for j in range(N)] if allowed(positions): print positions test = False
 * 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.

from mayavi.mlab import 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 code

=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). code format="python" import random, math, numpy N = 13; r = .25; R = 1./(1./r-1.); tmax = 1000; sig = .3; def dist(a,b): return math.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2) def allowed(positions): dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) print min(dists) return min(dists)>2*r def allowed_move(positions,k,newpos): dists = [dist(positions[l], newpos) for l in range(k)+range(k+1,N)] return min(dists)>2*r def unif_sphere: sigma = 1./math.sqrt(3.) x = [random.gauss(0,sigma) for i in range(3)] norm = math.sqrt(sum(xk**2. for xk in x)) return [xk/norm for xk in x] test = True while test: positions = [unif_sphere for j in range(N)] if allowed(positions): print positions test = False count = 0 for t in range(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 print "Rejection rate: ", 1.*count/tmax code
 * 1) initialisation
 * 1) moves

=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. code format="python" import random, math, numpy N = 13;r = .05;R = 1./(1./r-1.);tmax = 800000;sig = .25;gamma = 0.05 random.seed(137) def dist(a,b): return math.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2) def allowed(positions): dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) return min(dists)>2*r def allowed_move(positions,k,newpos): dists = [dist(positions[l], newpos) for l in range(k)+range(k+1,N)] return min(dists)>2*r def unif_sphere: sigma = 1./math.sqrt(3.) x = [random.gauss(0,sigma) for i in range(3)] norm = math.sqrt(sum(xk**2. for xk in x)) return [xk/norm for xk in x] test = True while test: positions = [unif_sphere for j in range(N)] if allowed(positions): test = False count = 0 for t in range(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 if t%100==0: if t%10000==0:print ' acceptance rate = ', 1.*count/(1+t) dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) Upsilon = min(dists)/2. r = r+gamma*(Upsilon-r)
 * 1) initialisation
 * 1) moves

R = 1./(1./r-1.) print 'density',1.*N/2.*(1.-math.sqrt(1-r*r)) print 'final r',r print 'final R',R

code

=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.

code format="python" import random, math, numpy N = 13; r = .2; R = 1./(1./r-1.); tmax = 500000; sig = .25; gamma = 0.1 random.seed(138) def dist(a,b): return math.sqrt((b[0]-a[0])**2+(b[1]-a[1])**2+(b[2]-a[2])**2) def allowed(positions): dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) return min(dists)>2*r def allowed_move(positions,k,newpos): dists = [dist(positions[l], newpos) for l in range(k)+range(k+1,N)] return min(dists)>2*r def unif_sphere: sigma = 1./math.sqrt(3.) x = [random.gauss(0,sigma) for i in range(3)] norm = math.sqrt(sum(xk**2. for xk in x)) return [xk/norm for xk in x] test = True while test: positions = [unif_sphere for j in range(N)] if allowed(positions): test = False count = 0 count2 = 0 for t in range(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+=1 if t%100 == 0: accept=count2/100. if accept<.2: sig = sig/2 if t%40000 == 0: print ' local acceptance rate = ', accept, '\tsigma = ', sig, '\tR = ', 1./(1./r-1.) count2 = 0 dists = [] for l in range(N): dists.extend([dist(positions[k],positions[l]) for k in range(l)]) Upsilon = min(dists)/2. r = r+gamma*(Upsilon-r)
 * 1) initialisation
 * 1) moves

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.mlab import 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 code =References=
 * K. Schütte und B. L. van der Waerden, //Das Problem der dreizehn Kugeln//, Mathematische Annalen 125, 325 (1953).
 * The kissing number problem.
 * G. Szpiro, Newton and the kissing problem.

url}?f=print| [Print this page]