Python/maths/dual_number_automatic_differentiation.py
Chris O 200429fc47
Dual Number Automatic Differentiation (#8760)
* Added dual_number_automatic_differentiation.py

* updating DIRECTORY.md

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update maths/dual_number_automatic_differentiation.py

---------

Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
2023-05-25 08:04:42 +02:00

142 lines
4.5 KiB
Python

from math import factorial
"""
https://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers
https://blog.jliszka.org/2013/10/24/exact-numeric-nth-derivatives.html
Note this only works for basic functions, f(x) where the power of x is positive.
"""
class Dual:
def __init__(self, real, rank):
self.real = real
if isinstance(rank, int):
self.duals = [1] * rank
else:
self.duals = rank
def __repr__(self):
return (
f"{self.real}+"
f"{'+'.join(str(dual)+'E'+str(n+1)for n,dual in enumerate(self.duals))}"
)
def reduce(self):
cur = self.duals.copy()
while cur[-1] == 0:
cur.pop(-1)
return Dual(self.real, cur)
def __add__(self, other):
if not isinstance(other, Dual):
return Dual(self.real + other, self.duals)
s_dual = self.duals.copy()
o_dual = other.duals.copy()
if len(s_dual) > len(o_dual):
o_dual.extend([1] * (len(s_dual) - len(o_dual)))
elif len(s_dual) < len(o_dual):
s_dual.extend([1] * (len(o_dual) - len(s_dual)))
new_duals = []
for i in range(len(s_dual)):
new_duals.append(s_dual[i] + o_dual[i])
return Dual(self.real + other.real, new_duals)
__radd__ = __add__
def __sub__(self, other):
return self + other * -1
def __mul__(self, other):
if not isinstance(other, Dual):
new_duals = []
for i in self.duals:
new_duals.append(i * other)
return Dual(self.real * other, new_duals)
new_duals = [0] * (len(self.duals) + len(other.duals) + 1)
for i, item in enumerate(self.duals):
for j, jtem in enumerate(other.duals):
new_duals[i + j + 1] += item * jtem
for k in range(len(self.duals)):
new_duals[k] += self.duals[k] * other.real
for index in range(len(other.duals)):
new_duals[index] += other.duals[index] * self.real
return Dual(self.real * other.real, new_duals)
__rmul__ = __mul__
def __truediv__(self, other):
if not isinstance(other, Dual):
new_duals = []
for i in self.duals:
new_duals.append(i / other)
return Dual(self.real / other, new_duals)
raise ValueError()
def __floordiv__(self, other):
if not isinstance(other, Dual):
new_duals = []
for i in self.duals:
new_duals.append(i // other)
return Dual(self.real // other, new_duals)
raise ValueError()
def __pow__(self, n):
if n < 0 or isinstance(n, float):
raise ValueError("power must be a positive integer")
if n == 0:
return 1
if n == 1:
return self
x = self
for _ in range(n - 1):
x *= self
return x
def differentiate(func, position, order):
"""
>>> differentiate(lambda x: x**2, 2, 2)
2
>>> differentiate(lambda x: x**2 * x**4, 9, 2)
196830
>>> differentiate(lambda y: 0.5 * (y + 3) ** 6, 3.5, 4)
7605.0
>>> differentiate(lambda y: y ** 2, 4, 3)
0
>>> differentiate(8, 8, 8)
Traceback (most recent call last):
...
ValueError: differentiate() requires a function as input for func
>>> differentiate(lambda x: x **2, "", 1)
Traceback (most recent call last):
...
ValueError: differentiate() requires a float as input for position
>>> differentiate(lambda x: x**2, 3, "")
Traceback (most recent call last):
...
ValueError: differentiate() requires an int as input for order
"""
if not callable(func):
raise ValueError("differentiate() requires a function as input for func")
if not isinstance(position, (float, int)):
raise ValueError("differentiate() requires a float as input for position")
if not isinstance(order, int):
raise ValueError("differentiate() requires an int as input for order")
d = Dual(position, 1)
result = func(d)
if order == 0:
return result.real
return result.duals[order - 1] * factorial(order)
if __name__ == "__main__":
import doctest
doctest.testmod()
def f(y):
return y**2 * y**4
print(differentiate(f, 9, 2))