Symbolic algebra for interaction potentials

[ajkluber]

I want to propose using sympy (a computer algebra system) for the interaction potentials, for the following reasons,

1. Using symbolic functions in sympy we can easily differentiate with respect to any variable or parameter. The symbolic functions can be transformed into regular python functions that return a number (see the example code below). This can straightforwardly replace some of the code in pairwise.py and remove the derivatives we did by hand (that are nasty). Not having to take derivatives by hand will also make development easier in the future.

import numpy as np
import sympy

# define symbolic variables
eps, r, r0, width = sympy.symbols("eps,r,r0,width")

# define symbolic function
G = eps*sympy.exp(-sympy.Rational(1, 2)*((r - r0)/width)**2)

# set fixed parameters
Gsub = G.subs([(eps, 1), (r0, 1), (width, 1)])

# Can evaluate the function at a vector of values
G_r = sympy.lambdify(r, Gsub, modules="numpy")
y = G_r(np.linspace(-1, 1, 1000))

# sympy allows for easy differentiation with respect to any variable
dGdr = sympy.diff(Gsub, r)
dGdr_r = sympy.lambdify(r, dGdr_sub, modules="numpy")

1. Along those same lines this fits with my long-term vision for how model_builder would interface with parameter fitting packages. We have moved ahead doing what has been the most logical and convenient at the time, but I think we have drifted from what I originally intended. My vision is that we should be able to fit any parameter that enters the Hamiltonian. Currently we are biased towards thinking about epsilons because it has been the easiest to think about, and this is completely natural and acceptable, but I want us to move in a more general direction. For example, we should be able to fit any parameter we want, regardless of where it enters the terms of the Hamiltonian (e.g. non-linear parameters). This is made much easier with sympy because of the symbolic differentiation support described above.

Maybe make this a milestone?