'Toolkit for automatic-differentiation of Python functions'
# Resources for automatic-differentiation:
# https://en.wikipedia.org/wiki/Automatic_differentiation
# https://justindomke.wordpress.com/2009/02/17/
# http://www.autodiff.org/
from __future__ import division
import math
## Python 2 vs 3 Compatibility Functions #################################
if not hasattr(math, 'isclose'):
math.isclose = lambda x, y, rel_tol=1e-09: abs(x/y - 1.0) <= rel_tol
if not hasattr(math, 'log2'):
math.log2 = lambda x: math.log(x) / math.log(2.0)
## Dual Number Class #####################################################
class Num(float):
''' The auto-differentiation number class works likes a float
for a function input, but all operations on that number
will concurrently compute the derivative.
Creating Nums
-------------
New numbers are created with: Num(x, dx)
Make constants (not varying with respect to x) with: Num(3.5)
Make variables (that vary with respect to x) with: Num(3.5, 1.0)
The short-cut for Num(3.5, 1.0) is: Var(3.5)
Accessing Nums
--------------
Convert a num back to a float with: float(n)
The derivative is accessed with: n.dx
Or with a the short-cut function: d(n)
Functions of One Variable
-------------------------
>>> f = lambda x: cos(2.5 * x) ** 3
>>> y = f(Var(1.5)) # Evaluate at x=1.5
>>> y # f(1.5)
-0.552497105486732
>>> y.dx # f'(1.5)
2.88631746797551
Partial Derivatives and Gradients of Multi-variable Functions
-------------------------------------------------------------
The tool can also be used to compute gradients of multivariable
functions by making one of the inputs variable and the keeping
the remaining inputs constant:
>>> f = lambda x, y: x*y + sin(x)
>>> f(Num(2.5), Num(3.5)) # Evaluate at (2.5, 3.5)
9.348472144103956
>>> d(f(Var(2.5), Num(3.5))) # Partial with respect to x
2.6988563844530664
>>> d(f(Num(2.5), Var(3.5))) # Partial with respect to y
2.5
>>> gradient(f, (2.5, 3.5))
(2.6988563844530664, 2.5)
See: https://www.wolframalpha.com/input/?lk=3&i=grad(x*y+%2B+sin(x))
'''
# Tables of Derivatives:
# http://hyperphysics.phy-astr.gsu.edu/hbase/math/derfunc.html
# http://tutorial.math.lamar.edu/pdf/Common_Derivatives_Integrals.pdf
# http://www.nps.edu/Academics/Schools/GSEAS/Departments/Math/pdf_sources/BlueBook27.pdf
# https://www.wolframalpha.com/input/?lk=3&i=d%2Fdx(u(x)%5E(v(x)))
__slots__ = ['dx']
def __new__(cls, value, dx=0.0):
if isinstance(value, cls): return value
inst = float.__new__(cls, value)
inst.dx = dx
return inst
def __add__(u, v):
return Num(float(u) + float(v), d(u) + d(v))
def __sub__(u, v):
return Num(float(u) - float(v), d(u) - d(v))
def __mul__(u, v):
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u * v, u * dv + v * du)
def __truediv__(u, v):
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u / v, (v * du - u * dv) / v ** 2.0)
def __pow__(u, v):
u, v, du, dv = float(u), float(v), d(u), d(v)
# XXX This doesn't handle negative u values when dv == 0.0
return Num(u ** v,
(v * u ** (v - 1.0) * du if du else 0.0) +
(math.log(u) * u ** v * dv if dv else 0.0))
def __mod__(u, v):
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u % v, du - u // v * dv)
def __pos__(u):
return u
def __neg__(u):
return Num(-float(u), -d(u))
__radd__ = __add__
__rmul__ = __mul__
def __rsub__(self, other):
return -(self - other)
def __rtruediv__(self, other):
return Num(other) / self
def __rpow__(self, other):
return Num(other) ** self
def __rmod__(u, v):
return Num(v) % u
def __abs__(self):
return self if self >= 0.0 else -self
## Convenience Functions #################################################
Var = lambda x: Num(x, 1.0)
d = lambda x: getattr(x, 'dx', 0.0)
## Math Module Functions and Constants ###################################
sqrt = lambda u: Num(math.sqrt(u), d(u) / (2.0 * math.sqrt(u)))
log = lambda u: Num(math.log(u), d(u) / float(u))
log2 = lambda u: Num(math.log2(u), d(u) / (float(u) * math.log(2.0)))
log10 = lambda u: Num(math.log10(u), d(u) / (float(u) * math.log(10.0)))
log1p = lambda u: Num(math.log1p(u), d(u) / (float(u) + 1.0))
exp = lambda u: Num(math.exp(u), math.exp(u) * d(u))
expm1 = lambda u: Num(math.expm1(u), math.exp(u) * d(u))
sin = lambda u: Num(math.sin(u), math.cos(u) * d(u))
cos = lambda u: Num(math.cos(u), -math.sin(u) * d(u))
tan = lambda u: Num(math.tan(u), d(u) / math.cos(u) ** 2.0)
sinh = lambda u: Num(math.sinh(u), math.cosh(u) * d(u))
cosh = lambda u: Num(math.cosh(u), math.sinh(u) * d(u))
tanh = lambda u: Num(math.tanh(u), d(u) / math.cosh(u) ** 2.0)
asin = lambda u: Num(math.asin(u), d(u) / math.sqrt(1.0 - float(u) ** 2.0))
acos = lambda u: Num(math.acos(u), -d(u) / math.sqrt(1.0 - float(u) ** 2.0))
atan = lambda u: Num(math.atan(u), d(u) / (1.0 + float(u) ** 2.0))
asinh = lambda u: Num(math.asinh(u), d(u) / math.hypot(u, 1.0))
acosh = lambda u: Num(math.acosh(u), d(u) / math.sqrt(float(u) ** 2.0 - 1.0))
atanh = lambda u: Num(math.atanh(u), d(u) / (1.0 - float(u) ** 2.0))
radians = lambda u: Num(math.radians(u), math.radians(d(u)))
degrees = lambda u: Num(math.degrees(u), math.degrees(d(u)))
hypot = lambda u, v: Num(math.hypot(u, v), (u * d(u) + v * d(v)) / math.hypot(u, v))
fsum = lambda u: Num(math.fsum(map(float, u)), math.fsum(map(d, u)))
copysign = lambda u, v: Num(math.copysign(u, v),
d(u) if math.copysign(1.0, float(u) * float(v)) > 0.0 else -d(u))
ceil = lambda u: Num(math.ceil(u), 0.0)
floor = lambda u: Num(math.floor(u), 0.0)
trunc = lambda u: Num(math.trunc(u), 0.0)
pi = Num(math.pi)
e = Num(math.e)
fabs = abs
fmod = lambda u, v: u % v
## Vector Functions ######################################################
def partial(func, point, index):
''' Partial derivative at a given point
>>> func = lambda x, y: x*y + sin(x)
>>> point = (2.5, 3.5)
>>> partial(func, point, 0) # Partial with respect to x
2.6988563844530664
>>> partial(func, point, 1) # Partial with respect to y
2.5
'''
return d(func(*[Num(x, i==index) for i, x in enumerate(point)]))
def gradient(func, point):
''' Vector of the partial derivatives of a scalar field
>>> func = lambda x, y: x*y + sin(x)
>>> point = (2.5, 3.5)
>>> gradient(func, point)
(2.6988563844530664, 2.5)
See: https://www.wolframalpha.com/input/?lk=3&i=grad(x*y+%2B+sin(x))
'''
return tuple(partial(func, point, index) for index in range(len(point)))
def directional_derivative(func, point, direction):
''' The dot product of the gradient and a direction vector.
Computed directly with a single function call.
>>> func = lambda x, y: x*y + sin(x)
>>> point = (2.5, 3.5)
>>> direction = (1.5, -2.2)
>>> directional_derivative(func, point, direction)
-1.4517154233204006
Same result as separately computing and dotting the gradient:
>>> math.fsum(g * d for g, d in zip(gradient(func, point), direction))
-1.4517154233204002
See: https://en.wikipedia.org/wiki/Directional_derivative
'''
return d(func(*map(Num, point, direction)))
def divergence(F, point):
''' Sum of the partial derivatives of a vector field
>>> F = lambda x, y, z: (x*y+sin(x)+3*x, x-y-5*x, cos(2*x)-sin(y)**2)
>>> divergence(F, (3.5, 2.1, -3.3))
3.163543312709203
# http://www.wolframalpha.com/input/?i=div+%7Bx*y%2Bsin(x)%2B3*x,+x-y-5*x,+cos(2*x)-sin(y)%5E2%7D
>>> x, y, z = (3.5, 2.1, -3.3)
>>> math.cos(x) + y + 2
3.1635433127092036
>>> F = lambda x, y, z: (8 * exp(-x), cosh(z), - y**2)
>>> divergence(F, (2, -1, 4))
-1.0826822658929016
# https://www.youtube.com/watch?v=S2rT2zK2bdo
>>> x, y, z = (2, -1, 4)
>>> -8 * math.exp(-x)
-1.0826822658929016
'''
return math.fsum(d(F(*[Num(x, i==index) for i, x in enumerate(point)])[index])
for index in range(len(point)))
def curl(F, point):
''' Rotation around a vector field
>>> F = lambda x, y, z: (x*y+sin(x)+3*x, x-y-5*x, cos(2*x)-sin(y)**2)
>>> curl(F, (3.5, 2.1, -3.3))
(0.8715757724135881, 1.3139731974375781, -7.5)
# http://www.wolframalpha.com/input/?i=curl+%7Bx*y%2Bsin(x)%2B3*x,+x-y-5*x,+cos(2*x)-sin(y)%5E2%7D
>>> x, y, z = (3.5, 2.1, -3.3)
>>> (-2 * math.sin(y) * math.cos(y), 2 * math.sin(2 * x), -x - 4)
(0.8715757724135881, 1.3139731974375781, -7.5)
# https://www.youtube.com/watch?v=UW4SQz29TDc
>>> F = lambda x, y, z: (y**4 - x**2 * z**2, x**2 + y**2, -x**2 * y * z)
>>> curl(F, (1, 3, -2))
(2.0, -8.0, -106.0)
>>> F = lambda x, y, z: (8 * exp(-x), cosh(z), - y**2)
>>> curl(F, (2, -1, 4))
(-25.289917197127753, 0.0, 0.0)
# https://www.youtube.com/watch?v=S2rT2zK2bdo
>>> x, y, z = (2, -1, 4)
>>> (-(x * y + math.sinh(z)), 0.0, 0.0)
(-25.289917197127753, 0.0, 0.0)
'''
x, y, z = point
_, Fyx, Fzx = map(d, F(Var(x), y, z))
Fxy, _, Fzy = map(d, F(x, Var(y), z))
Fxz, Fyz, _ = map(d, F(x, y, Var(z)))
return (Fzy - Fyz, Fxz - Fzx, Fyx - Fxy)
if __name__ == '__main__':
# River flow example: https://www.youtube.com/watch?v=vvzTEbp9lrc
W = 20 # width of river in meters
C = 0.1 # max flow divided by (W/2)**2
F = lambda x, y=0, z=0: (0.0, C * x * (W - x), 0.0)
for x in range(W+1):
print('%d --> %r' % (x, curl(F, (x, 0.0, 0.0))))
import doctest
import random
def numeric_derivative(func, x, eps=0.001):
'Estimate the derivative using numerical methods'
y0 = func(x - eps)
y1 = func(x + eps)
return (y1 - y0) / (2.0 * eps)
def test(x_array, *testcases):
for f in testcases:
print(f.__name__.center(40))
print('-' * 40)
for x in map(Var, x_array):
y = f(x)
actual = d(y)
expected = numeric_derivative(f, x, 2**-16)
print('%7.3f %12.4f %12.4f' % (x, actual, expected))
assert math.isclose(expected, actual, rel_tol=1e-5)
print('')
def test_pow_const_base(x):
return 3.1 ** (2.3 * x + 0.4)
def test_pow_const_exp(x):
return (2.3 * x + 0.4) ** (-1.3)
def test_pow_const_negative_base(x):
return (-3.1) ** (2.3 * x + 0.4)
def test_pow_general(x):
return (x / 3.5) ** sin(3.5 * x)
def test_hyperbolics(x):
return 3 * cosh(1/x) + 5 * sinh(x / 2.5) ** 2 - 0.7 * tanh(1.7 / x) ** 1.5
def test_sqrt(x):
return cos(sqrt(abs(sin(x) + 5)))
def test_conversions(x):
return degrees(x ** 1.5 + 18) * radians(0.83 ** x + 37)
def test_hypot(x):
return hypot(sin(x), cos(1.1 / x))
def test_rounders(x):
return tan(x) * floor(cos(x**2 + 0.37)*2.7) + \
log(x) * ceil(cos(x**3 + 0.31)*12.1) * 10.1 + \
exp(x) * trunc(sin(x**1.4 + 8.0)) * 1234.567
def test_inv_trig(x):
return atan((x - 0.303) ** 2.9 + 0.1234) + \
acos((x - 4.1) / 3.113) * 5 + \
asin((x - 4.3) / 3.717)
def test_mod(x):
return 137.1327 % (sin(x + .3) * 40.123) + cos(x) % 5.753
def test_logs(x):
return log2(fabs(sin(x))) + log10(fabs(cos(x))) + log1p(fabs(tan(x)))
def test_fsum(x):
n = 100
random.seed(8675309)
data = [Num(random.random()**x, random.random()**x) for i in range(n)]
return fsum(data)
def test_inv_hyperbolics(x):
return (acosh(x ** 1.1234567 + 0.89) + 3.51 * asinh(x ** 1.234567 + 8.9) +
atanh(x / 15.0823))
def test_copysign(x):
return (copysign(7.17 * x + 5.11, 1.0) + copysign(4.1 * x, 0.0) +
copysign(8.909 * x + 0.18, -0.0) + copysign(4.321 * x + .12, -1.0) +
copysign(-3.53 * x + 11.5, 1.0) + copysign(-1.4 * x + 2.1, 0.0) +
copysign(-9.089 * x + 0.813, -0.0) + copysign(-1.2347 * x, -1.0) +
copysign(sin(x), x - math.pi))
def test_combined(x):
return 1.7 - 3 * cos(x) ** 2 / sin(3 * x) * 0.1 * exp(+cos(x)) + \
sqrt(abs(x - 4.13)) + tan(2.5 * x) * log(3.1 * x**1.5) + \
(4.7 * x + 3.1) ** cos(.43*x + 8.1) - 2.9 + tan(-x) + \
sqrt(radians(log(x) + 1.7)) + e / x + expm1(x / pi)
x_array = [2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]
working = [test_combined, test_pow_const_base, test_pow_const_exp,
test_pow_general, test_hyperbolics, test_sqrt, test_copysign,
test_inv_trig, test_conversions, test_hypot, test_rounders,
test_inv_hyperbolics, test_mod, test_logs, test_fsum]
broken = [test_pow_const_negative_base]
test(x_array, *working)
print(doctest.testmod())
Diff to Previous Revision
--- revision 2 2016-02-19 05:15:00
+++ revision 3 2016-02-21 21:49:56
@@ -80,25 +80,20 @@
return inst
def __add__(u, v):
- v = Num(v)
return Num(float(u) + float(v), d(u) + d(v))
def __sub__(u, v):
- v = Num(v)
return Num(float(u) - float(v), d(u) - d(v))
def __mul__(u, v):
- v = Num(v)
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u * v, u * dv + v * du)
def __truediv__(u, v):
- v = Num(v)
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u / v, (v * du - u * dv) / v ** 2.0)
def __pow__(u, v):
- v = Num(v)
u, v, du, dv = float(u), float(v), d(u), d(v)
# XXX This doesn't handle negative u values when dv == 0.0
return Num(u ** v,
@@ -106,7 +101,6 @@
(math.log(u) * u ** v * dv if dv else 0.0))
def __mod__(u, v):
- v = Num(v)
u, v, du, dv = float(u), float(v), d(u), d(v)
return Num(u % v, du - u // v * dv)
@@ -136,7 +130,7 @@
## Convenience Functions #################################################
Var = lambda x: Num(x, 1.0)
-d = lambda x: x.dx
+d = lambda x: getattr(x, 'dx', 0.0)
## Math Module Functions and Constants ###################################
sqrt = lambda u: Num(math.sqrt(u), d(u) / (2.0 * math.sqrt(u)))
@@ -167,8 +161,8 @@
ceil = lambda u: Num(math.ceil(u), 0.0)
floor = lambda u: Num(math.floor(u), 0.0)
trunc = lambda u: Num(math.trunc(u), 0.0)
-pi = Num(math.pi, 0.0)
-e = Num(math.e, 0.0)
+pi = Num(math.pi)
+e = Num(math.e)
fabs = abs
fmod = lambda u, v: u % v
@@ -272,9 +266,9 @@
'''
x, y, z = point
- _, Fyx, Fzx = map(d, F(Var(x), Num(y), Num(z)))
- Fxy, _, Fzy = map(d, F(Num(x), Var(y), Num(z)))
- Fxz, Fyz, _ = map(d, F(Num(x), Num(y), Var(z)))
+ _, Fyx, Fzx = map(d, F(Var(x), y, z))
+ Fxy, _, Fzy = map(d, F(x, Var(y), z))
+ Fxz, Fyz, _ = map(d, F(x, y, Var(z)))
return (Fzy - Fyz, Fxz - Fzx, Fyx - Fxy)
@@ -282,7 +276,101 @@
# River flow example: https://www.youtube.com/watch?v=vvzTEbp9lrc
W = 20 # width of river in meters
- C = .1 # max flow divided by (W/2)**2
- F = lambda x, y=0, z=0: (Num(0), C * x * (W - x), Num(0))
+ C = 0.1 # max flow divided by (W/2)**2
+ F = lambda x, y=0, z=0: (0.0, C * x * (W - x), 0.0)
for x in range(W+1):
print('%d --> %r' % (x, curl(F, (x, 0.0, 0.0))))
+
+ import doctest
+ import random
+
+ def numeric_derivative(func, x, eps=0.001):
+ 'Estimate the derivative using numerical methods'
+ y0 = func(x - eps)
+ y1 = func(x + eps)
+ return (y1 - y0) / (2.0 * eps)
+
+ def test(x_array, *testcases):
+ for f in testcases:
+ print(f.__name__.center(40))
+ print('-' * 40)
+ for x in map(Var, x_array):
+ y = f(x)
+ actual = d(y)
+ expected = numeric_derivative(f, x, 2**-16)
+ print('%7.3f %12.4f %12.4f' % (x, actual, expected))
+ assert math.isclose(expected, actual, rel_tol=1e-5)
+ print('')
+
+ def test_pow_const_base(x):
+ return 3.1 ** (2.3 * x + 0.4)
+
+ def test_pow_const_exp(x):
+ return (2.3 * x + 0.4) ** (-1.3)
+
+ def test_pow_const_negative_base(x):
+ return (-3.1) ** (2.3 * x + 0.4)
+
+ def test_pow_general(x):
+ return (x / 3.5) ** sin(3.5 * x)
+
+ def test_hyperbolics(x):
+ return 3 * cosh(1/x) + 5 * sinh(x / 2.5) ** 2 - 0.7 * tanh(1.7 / x) ** 1.5
+
+ def test_sqrt(x):
+ return cos(sqrt(abs(sin(x) + 5)))
+
+ def test_conversions(x):
+ return degrees(x ** 1.5 + 18) * radians(0.83 ** x + 37)
+
+ def test_hypot(x):
+ return hypot(sin(x), cos(1.1 / x))
+
+ def test_rounders(x):
+ return tan(x) * floor(cos(x**2 + 0.37)*2.7) + \
+ log(x) * ceil(cos(x**3 + 0.31)*12.1) * 10.1 + \
+ exp(x) * trunc(sin(x**1.4 + 8.0)) * 1234.567
+
+ def test_inv_trig(x):
+ return atan((x - 0.303) ** 2.9 + 0.1234) + \
+ acos((x - 4.1) / 3.113) * 5 + \
+ asin((x - 4.3) / 3.717)
+
+ def test_mod(x):
+ return 137.1327 % (sin(x + .3) * 40.123) + cos(x) % 5.753
+
+ def test_logs(x):
+ return log2(fabs(sin(x))) + log10(fabs(cos(x))) + log1p(fabs(tan(x)))
+
+ def test_fsum(x):
+ n = 100
+ random.seed(8675309)
+ data = [Num(random.random()**x, random.random()**x) for i in range(n)]
+ return fsum(data)
+
+ def test_inv_hyperbolics(x):
+ return (acosh(x ** 1.1234567 + 0.89) + 3.51 * asinh(x ** 1.234567 + 8.9) +
+ atanh(x / 15.0823))
+
+ def test_copysign(x):
+ return (copysign(7.17 * x + 5.11, 1.0) + copysign(4.1 * x, 0.0) +
+ copysign(8.909 * x + 0.18, -0.0) + copysign(4.321 * x + .12, -1.0) +
+ copysign(-3.53 * x + 11.5, 1.0) + copysign(-1.4 * x + 2.1, 0.0) +
+ copysign(-9.089 * x + 0.813, -0.0) + copysign(-1.2347 * x, -1.0) +
+ copysign(sin(x), x - math.pi))
+
+ def test_combined(x):
+ return 1.7 - 3 * cos(x) ** 2 / sin(3 * x) * 0.1 * exp(+cos(x)) + \
+ sqrt(abs(x - 4.13)) + tan(2.5 * x) * log(3.1 * x**1.5) + \
+ (4.7 * x + 3.1) ** cos(.43*x + 8.1) - 2.9 + tan(-x) + \
+ sqrt(radians(log(x) + 1.7)) + e / x + expm1(x / pi)
+
+ x_array = [2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]
+ working = [test_combined, test_pow_const_base, test_pow_const_exp,
+ test_pow_general, test_hyperbolics, test_sqrt, test_copysign,
+ test_inv_trig, test_conversions, test_hypot, test_rounders,
+ test_inv_hyperbolics, test_mod, test_logs, test_fsum]
+ broken = [test_pow_const_negative_base]
+ test(x_array, *working)
+
+ print(doctest.testmod())
