mirror of
https://github.com/TheAlgorithms/Python.git
synced 2024-10-05 21:29:29 +00:00
4b79d771cd
* Add more ruff rules * Add more ruff rules * pre-commit: Update ruff v0.0.269 -> v0.0.270 * Apply suggestions from code review * Fix doctest * Fix doctest (ignore whitespace) * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: Dhruv Manilawala <dhruvmanila@gmail.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
142 lines
4.5 KiB
Python
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))
|