public
anonymous / Steepest Descent
Created

Algorithm to place N points on a sphere with equal distances to all their closest neighbours. Uses the golden spiral algorithm to place the initial points followed by energy minimization using the steepest descent method to optimize the distances.

  • Download Gist
Steepest Descent
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
#!/usr/bin/env python
 
"""
Generates a set of points on a sphere and evenly distributes them
by minimizing the potential energy according to a 'Coloumb' potential.
 
Joao Rodrigues, 2013
"""
 
import numpy
from scipy.spatial import distance
 
def generate_points(N):
"""
Generates a distribution of N points on a sphere based
on the Golden Spiral Algorithm.
"""
inc = numpy.pi * (3 - numpy.sqrt(5))
off = 2.0 / N
pts = [(numpy.cos(k * inc)*numpy.sqrt(1.0 - (k * off - 1 + (off / 2.0))**2), (k * off - 1 + (off / 2.0)), numpy.sin(k * inc)*numpy.sqrt(1.0 - (k * off - 1 + (off / 2.0))**2)) for k in xrange(N)]
return numpy.array(pts, dtype=numpy.float64)
 
def calc_coulomb(coords):
"""
Ecoulomb = E(1/r)
"""
distances = distance.pdist(coords, 'euclidean')
return numpy.power(distances, -1).sum()
def calc_forces(system):
"""
Calculates forces acting on all particles.
"""
forces = [ [0,0,0] for particle in system ]
nparticles = len(system)
for i in xrange(nparticles):
particle_i = system[i]
for j in xrange(i+1, nparticles):
particle_j = system[j]
r = particle_i - particle_j
f = numpy.linalg.norm(r)**-3
forces[i][0] += r[0]*f
forces[i][1] += r[1]*f
forces[i][2] += r[2]*f
 
forces[j][0] -= r[0]*f
forces[j][1] -= r[1]*f
forces[j][2] -= r[2]*f
return numpy.array(forces, dtype=numpy.float64)
 
def steepest_descent(coordinates, max_steps, step_size, min_step_size, min_energy_delta, verbose=True):
"""
Steepest Descent minimization algorithm
 
http://mathworld.wolfram.com/MethodofSteepestDescent.html
"""
 
nparticles = len(coordinates)
for k in xrange(options.nsteps):
# Backup coordinates
old_coordinates = coordinates.copy()
# 3.1 Calculate energy
cur_energy = calc_coulomb(old_coordinates)
 
# 3.2 Calculate forces
forces = calc_forces(coordinates)
force_normal = numpy.linalg.norm(forces)
 
# 3.3 Move particles
if force_normal > 0:
for p in xrange(nparticles):
directions = (forces[p]/force_normal)*step_size
coordinates[p] += directions
coordinates[p] /= numpy.linalg.norm(coordinates[p])
 
# 3.4 Check new energy and accept move if favourable + accelerate in that direction
new_energy = calc_coulomb(coordinates)
if new_energy >= cur_energy:
step_size /= 2
if step_size < min_step_size:
print "-- Step Size Converged (%1.5f < %1.5f)" %(step_size, min_step_size)
break
coordinates = old_coordinates
else:
if verbose:
print '-- Step %i: Energy %12.8f' %(k+1, cur_energy)
delta_e = abs(new_energy - cur_energy)
if delta_e < min_energy_delta:
print "-- Energy Converged (%12.12f < %12.12f)" %(delta_e, min_energy_delta)
break
step_size *= 2
cur_energy = new_energy
return (coordinates, new_energy)
 
def plot_coordinates(initial, final):
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x0, y0, z0 = zip(*initial)
ax.scatter(x0, y0, z0, zdir='z', s=20, c='b')
xmin, ymin, zmin = zip(*final)
ax.scatter(xmin, ymin, zmin, zdir='z', s=20, c='r')
u, v = numpy.mgrid[0:2*numpy.pi:20j, 0:numpy.pi:10j]
x=numpy.cos(u)*numpy.sin(v)
y=numpy.sin(u)*numpy.sin(v)
z=numpy.cos(v)
ax.plot_wireframe(x, y, z, color="gray")
plt.show()
if __name__ == "__main__":
import sys
 
from optparse import OptionParser
parser = OptionParser(usage="%prog <number of points> [options] [--help]")
parser.add_option('-d', action="store", dest="step_size", type="float",
default=0.01, help="Step size for the minimizer")
parser.add_option('-n', action="store", dest="nsteps", type="int",
default=1000, help="Number of steps to minimize")
parser.add_option('-c', action="store", dest="n_convergence", type="float",
default=(10**-5), help="Minimum Step Size (Convergence Criteria)")
parser.add_option('-e', action="store", dest="e_convergence", type="float",
default=(10**-10), help="Minimum Energy Change (Convergence Criteria)")
parser.add_option('-q', action="store_true", dest="quiet",
help="Supresses output from minimization algorithm")
parser.add_option('-p', action="store_true", dest="plot",
help="Generates a plot of the initial and final coordinates [req. matplotlib]")
options, args = parser.parse_args()
 
# 0. Input
try:
npoints = int(args[0])
except (IndexError, ValueError):
parser.print_usage()
sys.exit(1)
# 1. Generate Spiral Points
initial_coordinates = generate_points(npoints)
# 2. Minimization (steepest descent)
optimized_coordinates, final_energy = steepest_descent(initial_coordinates.copy(),
options.nsteps,
options.step_size,
options.n_convergence,
options.e_convergence,
verbose=(not options.quiet))
 
print '-- Final Energy: %12.8f' %(final_energy)
 
# Plot initial vs final positions on the sphere
if options.plot:
plot_coordinates(initial_coordinates, optimized_coordinates)
 
sys.exit(0)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.