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 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 [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)
to join this conversation on GitHub. Already have an account? Sign in to comment