Create a gist now

Instantly share code, notes, and snippets.

anonymous /Steepest Descent
Created Feb 2, 2013

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.
#!/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
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)
coordinates = old_coordinates
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)
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]
ax.plot_wireframe(x, y, z, color="gray")
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
npoints = int(args[0])
except (IndexError, ValueError):
# 1. Generate Spiral Points
initial_coordinates = generate_points(npoints)
# 2. Minimization (steepest descent)
optimized_coordinates, final_energy = steepest_descent(initial_coordinates.copy(),
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment