dual_number_automatic_differentiation.py

frommathimportfactorial

"""
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.
"""


classDual:
def__init__(self, real, rank):
self.real=real
ifisinstance(rank, int):
self.duals= [1] *rank
else:
self.duals=rank

def__repr__(self):
s="+".join(f"{dual}E{n}"forn, dualinenumerate(self.duals, 1))
returnf"{self.real}+{s}"

defreduce(self):
cur=self.duals.copy()
whilecur[-1] ==0:
cur.pop(-1)
returnDual(self.real, cur)

def__add__(self, other):
ifnotisinstance(other, Dual):
returnDual(self.real+other, self.duals)
s_dual=self.duals.copy()
o_dual=other.duals.copy()
iflen(s_dual) >len(o_dual):
o_dual.extend([1] * (len(s_dual) -len(o_dual)))
eliflen(s_dual) <len(o_dual):
s_dual.extend([1] * (len(o_dual) -len(s_dual)))
new_duals= []
foriinrange(len(s_dual)):
new_duals.append(s_dual[i] +o_dual[i])
returnDual(self.real+other.real, new_duals)

__radd__=__add__

def__sub__(self, other):
returnself+other*-1

def__mul__(self, other):
ifnotisinstance(other, Dual):
new_duals= []
foriinself.duals:
new_duals.append(i*other)
returnDual(self.real*other, new_duals)
new_duals= [0] * (len(self.duals) +len(other.duals) +1)
fori, iteminenumerate(self.duals):
forj, jteminenumerate(other.duals):
new_duals[i+j+1] +=item*jtem
forkinrange(len(self.duals)):
new_duals[k] +=self.duals[k] *other.real
forindexinrange(len(other.duals)):
new_duals[index] +=other.duals[index] *self.real
returnDual(self.real*other.real, new_duals)

__rmul__=__mul__

def__truediv__(self, other):
ifnotisinstance(other, Dual):
new_duals= []
foriinself.duals:
new_duals.append(i/other)
returnDual(self.real/other, new_duals)
raiseValueError

def__floordiv__(self, other):
ifnotisinstance(other, Dual):
new_duals= []
foriinself.duals:
new_duals.append(i//other)
returnDual(self.real//other, new_duals)
raiseValueError

def__pow__(self, n):
ifn<0orisinstance(n, float):
raiseValueError("power must be a positive integer")
ifn==0:
return1
ifn==1:
returnself
x=self
for_inrange(n-1):
x*=self
returnx


defdifferentiate(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
 """
ifnotcallable(func):
raiseValueError("differentiate() requires a function as input for func")
ifnotisinstance(position, (float, int)):
raiseValueError("differentiate() requires a float as input for position")
ifnotisinstance(order, int):
raiseValueError("differentiate() requires an int as input for order")
d=Dual(position, 1)
result=func(d)
iforder==0:
returnresult.real
returnresult.duals[order-1] *factorial(order)


if__name__=="__main__":
importdoctest

doctest.testmod()

deff(y):
returny**2*y**4

print(differentiate(f, 9, 2))