Dual Number Automatic Differentiation

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))

DIRECTORY.md

@@ -549,6 +549,7 @@
  * [Dodecahedron](maths/dodecahedron.py)
  * [Double Factorial Iterative](maths/double_factorial_iterative.py)
  * [Double Factorial Recursive](maths/double_factorial_recursive.py)
+ * [Dual Number Automatic Differentiation](maths/dual_number_automatic_differentiation.py)
  * [Entropy](maths/entropy.py)
  * [Euclidean Distance](maths/euclidean_distance.py)
  * [Euclidean Gcd](maths/euclidean_gcd.py)

maths/dual_number_automatic_differentiation.py

@@ -0,0 +1,141 @@
+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))