add exploration for precision of the langevin function 1/tanh(x) + 1/x

explore/precision.py

@@ -190,13 +190,13 @@ def plotdiff(x, target, actual, label, diff):
     if diff == "relative":
         err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     elif diff == "absolute":
         err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     else:
         limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
+        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label, alpha=0.7)
 
 def make_ocl(function, name, source=[]):
     class Kernel(object):
@@ -412,6 +412,54 @@ def add_function(name, mp_function, np_function, ocl_function,
     np_function=lambda x: np.fmod(x, 2*np.pi),
     ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
 )
+
+def sas_langevin(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    #cutoff *= 10
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    # 4 terms gets to 1e-7 single, 1e-14 double. Can get to 1e-15 double by adding
+    # another 4 terms and setting cutoff at 1.0. Not worthwhile. Instead we would
+    # need an expansion about x somewhere between 1 and 10 for the interval [0.1, 100.]
+    #f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9. + xpsq/(11.0 + xpsq/(13. + xpsq/(15. + xpsq/17.)))))))
+    xp = x[~index]
+    f[~index] = 1/np.tanh(xp) - 1/xp
+    return f[0] if scalar else f
+
+def sas_langevin_x(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = 1. / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    xp = x[~index]
+    f[~index] = (1/np.tanh(xp) - 1/xp)/xp
+    return f[0] if scalar else f
+
+add_function(
+    name="langevin(x)",
+    mp_function=lambda x: (1/mp.tanh(x) - 1/x),
+    np_function=sas_langevin,
+    #ocl_function=make_ocl("return q < 0.7 ? q*(1./3. + q*q*(-1./45. + q*q*(2./945. + q*q*(-1./4725.) + q*q*(2./93555.)))) : 1/tanh(q) - 1/q;", "sas_langevin"),
+    ocl_function=make_ocl("return q < 1e-5 ? q/3. : 1/tanh(q) - 1/q;", "sas_langevin"),
+)
+add_function(
+    name="langevin(x)/x",
+    mp_function=lambda x: (1/mp.tanh(x) - 1/x)/x,
+    #np_function=lambda x: sas_langevin(x)/x, # Note: need to test for x=0
+    np_function=sas_langevin_x,
+    ocl_function=make_ocl("return q < 1e-5 ? 1./3. : (1/tanh(q) - 1/q)/q;", "sas_langevin_x"),
+)
 add_function(
     name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,