Welcome to ableiter
This is the documentation of the
package ableiter: A tool for the
numerical differentiation of functions
with several variables using
finite differences.
Installation
ableiter can be installed by
typing the following command
from within the main directory:
pip install .
Issues
In case you notice any bug or otherwise undesired behaviour, please contact me or raise an issue in the gitlab repository.
User guide
The package is able to perform
partial derivatives up to order
four with respect to any of the
variables of the function.
First-order derivatives can be
computed with the help of the
class First, second-order
derivatives with the class
Second, and so on. The usage
of these classes is documented
below.
In order to use the package
ableiter, it is required
to define the function to
be derived in the following way.
The function (called f from hereon)
should have one single argument.
This argument shall be a numpy
array with one dimension. The index
of the array runs over the variables
of f, i.e. x[0] is the first variable
of f, x[1] is the second variable,
and so on. The length of the input
array has to be equal to the number
of variables of f.
For instance, to define the function
one should type:
def f(x):
return x[0]**3 - x[1]**2 + x[0]**3 * x[1]**4 + x[1] * np.log(x[0])
First
ableiter.first.First
This class can be used to compute the partial first-order derivatives of a multivariate function.
__init__(self, f, n)
special
Sets the function f as the function whose partial derivatives are calculated, and n is the dimension of the input array of f.
Source code in ableiter/first.py
def __init__(self, f, n):
"""
Sets the function f as the function whose
partial derivatives are calculated, and n
is the dimension of the input array of f.
"""
self.f = f
self.n = n
self.x_eps = 1e-8
self.grf = np.zeros(shape=(n, ))
df_dxi(self, y, i)
Return the value of the first derivative with respect to the i-th variable in the point y.
Source code in ableiter/first.py
def df_dxi(self, y, i):
"""
Return the value of the first derivative with
respect to the i-th variable in the point y.
"""
x = np.copy(y)
if len(x) != self.n:
raise ValueError
f = abs(self.f(x))
if f > 1e0:
x_eps = self.x_eps * f**0.5
else:
x_eps = self.x_eps
xiold = x[i]
x[i] = xiold - x_eps
f1 = self.f(x)
x[i] = xiold + x_eps
f2 = self.f(x)
return 0.5 * (f2 - f1) / x_eps
gradf(self, x)
Return the gradient in the point x as a numpy array with length n.
Source code in ableiter/first.py
def gradf(self, x):
"""
Return the gradient in the point x
as a numpy array with length n.
"""
for i in range(0, self.n):
self.grf[i] = self.df_dxi(x, i)
return self.grf
Second
ableiter.second.Second
This class can be used to compute partial second-order derivatives of a multivariate function.
__init__(self, f, n)
special
Sets the function f as the function whose partial derivatives are calculated, and n is the dimension of the input array of f.
Source code in ableiter/second.py
def __init__(self, f, n):
"""
Sets the function f as the function whose
partial derivatives are calculated, and n
is the dimension of the input array of f.
"""
self.f = f
self.n = n
self.x_eps = 1e-6
d2f_dxidxj(self, y, i, j)
Return the value of the second derivative in the point y, where first the function is derived with respect to the i-th variable, and then the functions is derived with respect to the j-th variable.
Source code in ableiter/second.py
def d2f_dxidxj(self, y, i, j):
"""
Return the value of the second derivative
in the point y,
where first the function is derived with
respect to the i-th variable, and then
the functions is derived with respect
to the j-th variable.
"""
x = np.copy(y)
if len(x) != self.n:
raise ValueError
f = abs(self.f(x))
if f > 1e0:
x_eps = self.x_eps * f**0.5
else:
x_eps = self.x_eps
if i != j:
xiold = x[i]
xjold = x[j]
x[i] = xiold - x_eps
x[j] = xjold - x_eps
fmm = self.f(x)
# x[i] = xiold - x_eps
x[j] = xjold + x_eps
fmp = self.f(x)
x[i] = xiold + x_eps
# x[j] = xjold + x_eps
fpp = self.f(x)
# x[i] = xiold + x_eps
x[j] = xjold - x_eps
fpm = self.f(x)
y = 0.25 * (fpp - fpm - fmp + fmm) / x_eps**2
else:
xiold = x[i]
f0 = self.f(x)
x[i] = xiold - x_eps
fm = self.f(x)
x[i] = xiold + x_eps
fp = self.f(x)
y = -(2. * f0 - fm - fp) / x_eps**2
return y
Third
ableiter.third.Third
This class can be used to compute partial third-order derivatives of a multivariate function.
__init__(self, f, n)
special
Sets the function f as the function whose partial derivatives are calculated, and n is the dimension of the input array of f.
Source code in ableiter/third.py
def __init__(self, f, n):
"""
Sets the function f as the function whose
partial derivatives are calculated, and n
is the dimension of the input array of f.
"""
self.f = f
self.n = n
self.x_eps = 1e-4
d3f_dxidxjdxk(self, y, i, j, k)
Return the value of the third derivative in the point y, where first the function is derived with respect to the i-th variable, then the functions is derived with respect to the j-th variable, and then the function is derived with respect to the k-th variable.
Source code in ableiter/third.py
def d3f_dxidxjdxk(self, y, i, j, k):
"""
Return the value of the third derivative
in the point y,
where first the function is derived with
respect to the i-th variable, then
the functions is derived with respect
to the j-th variable, and then the
function is derived with respect to
the k-th variable.
"""
x = np.copy(y)
if len(x) != self.n:
raise ValueError
f = abs(self.f(x))
if f > 1e0:
x_eps = self.x_eps * f**0.5
else:
x_eps = self.x_eps
if (i == j) and (j == k):
xiold = x[i]
x[i] = xiold - 2. * x_eps
fmm = self.f(x)
x[i] = xiold - x_eps
fm0 = self.f(x)
x[i] = xiold + x_eps
fp0 = self.f(x)
x[i] = xiold + 2. * x_eps
fpp = self.f(x)
y = 0.5 * (2. * fm0 - 2. * fp0 - fmm + fpp) / \
x_eps**3
elif (i == j) or (j == k) or (i == k):
if i == j:
il = i
ir = k
elif j == k:
il = j
ir = i
else:
il = i
ir = j
xilold = x[il]
xirold = x[ir]
x[il] = xilold + 2. * x_eps
x[ir] = xirold + x_eps
f21 = self.f(x)
x[il] = xilold
x[ir] = xirold + x_eps
f01 = self.f(x)
x[il] = xilold - 2. * x_eps
x[ir] = xirold + x_eps
fm21 = self.f(x)
x[il] = xilold + 2. * x_eps
x[ir] = xirold - x_eps
f2m1 = self.f(x)
x[il] = xilold
x[ir] = xirold - x_eps
f0m1 = self.f(x)
x[il] = xilold - 2. * x_eps
x[ir] = xirold - x_eps
fm2m1 = self.f(x)
y = 0.125 * (f21 - 2. * f01 + fm21 - f2m1 + \
2. * f0m1 - fm2m1) / x_eps**3
else:
xiold = x[i]
xjold = x[j]
xkold = x[k]
x[i] = xiold + x_eps
x[j] = xjold + x_eps
x[k] = xkold + x_eps
fppp = self.f(x)
x[i] = xiold - x_eps
x[j] = xjold + x_eps
x[k] = xkold + x_eps
fmpp = self.f(x)
x[i] = xiold + x_eps
x[j] = xjold - x_eps
x[k] = xkold + x_eps
fpmp = self.f(x)
x[i] = xiold - x_eps
x[j] = xjold - x_eps
x[k] = xkold + x_eps
fmmp = self.f(x)
x[i] = xiold + x_eps
x[j] = xjold + x_eps
x[k] = xkold - x_eps
fppm = self.f(x)
x[i] = xiold - x_eps
x[j] = xjold + x_eps
x[k] = xkold - x_eps
fmpm = self.f(x)
x[i] = xiold + x_eps
x[j] = xjold - x_eps
x[k] = xkold - x_eps
fpmm = self.f(x)
x[i] = xiold - x_eps
x[j] = xjold - x_eps
x[k] = xkold - x_eps
fmmm = self.f(x)
y = 0.125 * (fppp - fmpp - fpmp + fmmp
- fppm + fmpm + fpmm - fmmm) / x_eps**3
return y
Fourth
ableiter.fourth.Fourth
This class can be used to compute partial fourth-order derivatives of a multivariate function.
__init__(self, f, n)
special
Sets the function f as the function whose partial derivatives are calculated, and n is the dimension of the input array of f.
Source code in ableiter/fourth.py
def __init__(self, f, n):
"""
Sets the function f as the function whose
partial derivatives are calculated, and n
is the dimension of the input array of f.
"""
self.f = f
self.n = n
self.x_eps = 1e-3
d4f_dxidxjdxkdxl(self, y, i, j, k, l)
Return the value of the fourth derivative in the point y, where first the function is derived with respect to the i-th variable, then the functions is derived with respect to the j-th variable, then the function is derived with respect to the k-th variable, and then the function is derived with respect to the l-th variable.
Source code in ableiter/fourth.py
def d4f_dxidxjdxkdxl(self, y, i, j, k, l):
"""
Return the value of the fourth derivative
in the point y,
where first the function is derived with
respect to the i-th variable, then
the functions is derived with respect
to the j-th variable, then the
function is derived with respect to
the k-th variable, and then the function
is derived with respect to the l-th variable.
"""
x = np.copy(y)
if len(x) != self.n:
raise ValueError
f = abs(self.f(x))
if f > 1e0:
x_eps = self.x_eps * f**0.5
else:
x_eps = self.x_eps
if len(set([i, j, k, l])) == 1:
complexity = 1
elif (i == j) and (k == l) and (i != k):
il = i
ir = k
complexity = 2
elif (i == k) and (j == l) and (i != j):
il = i
ir = j
complexity = 2
elif (i == l) and (j == k) and (i != j):
il = i
ir = j
complexity = 2
elif (i == j) and (i == k) and (i != l):
il = i
ir = l
complexity = 3
elif (i == j) and (i == l) and (i != k):
il = i
ir = k
complexity = 3
elif (i == k) and (i == l) and (i != j):
il = i
ir = j
complexity = 3
elif (i != j) and (j == k) and (j == l):
il = j
ir = i
complexity = 3
elif len(set([i, j, k, l])) == 4:
complexity = 4
elif (i == j) and (len(set([i, j, k, l])) == 3):
il = i
im = k
ir = l
complexity = 5
elif (i == k) and (len(set([i, j, k, l])) == 3):
il = i
im = j
ir = l
complexity = 5
elif (i == l) and (len(set([i, j, k, l])) == 3):
il = i
im = j
ir = k
complexity = 5
elif (j == k) and (len(set([i, j, k, l])) == 3):
il = j
im = i
ir = l
complexity = 5
elif (j == l) and (len(set([i, j, k, l])) == 3):
il = j
im = i
ir = k
complexity = 5
elif (k == l) and (len(set([i, j, k, l])) == 3):
il = k
im = i
ir = j
complexity = 5
else:
raise NotImplementedError
if complexity == 2:
xilold = x[il]
xirold = x[ir]
x[il] = xilold + x_eps
x[ir] = xirold + x_eps
f11 = self.f(x)
x[il] = xilold
x[ir] = xirold + x_eps
f01 = self.f(x)
x[il] = xilold - x_eps
x[ir] = xirold + x_eps
fm11 = self.f(x)
x[il] = xilold + x_eps
x[ir] = xirold
f10 = self.f(x)
x[il] = xilold
x[ir] = xirold
f00 = self.f(x)
x[il] = xilold - x_eps
x[ir] = xirold
fm10 = self.f(x)
x[il] = xilold + x_eps
x[ir] = xirold - x_eps
f1m1 = self.f(x)
x[il] = xilold
x[ir] = xirold - x_eps
f0m1 = self.f(x)
x[il] = xilold - x_eps
x[ir] = xirold - x_eps
fm1m1 = self.f(x)
y = (f11 - 2. * f01 + fm11
- 2. * f10 + 4. * f00 - 2. * fm10
+ f1m1 - 2. * f0m1 + fm1m1) / \
x_eps**4
elif complexity == 1:
xiold = x[i]
f0 = self.f(x)
x[i] = xiold + 2. * x_eps
f2 = self.f(x)
x[i] = xiold + x_eps
f1 = self.f(x)
x[i] = xiold - x_eps
fm1 = self.f(x)
x[i] = xiold - 2. * x_eps
fm2 = self.f(x)
y = (f2 - 4. * f1 + 6. * f0
- 4. * fm1 + fm2) / x_eps**4
elif complexity == 3:
xilold = x[il]
xirold = x[ir]
x[il] = xilold + 3. * x_eps
x[ir] = xirold + x_eps
f13 = self.f(x)
# x[il] = xilold + 3. * x_eps
x[ir] = xirold - x_eps
fm13 = self.f(x)
x[il] = xilold + x_eps
x[ir] = xirold + x_eps
f11 = self.f(x)
# x[il] = xilold + x_eps
x[ir] = xirold - x_eps
fm11 = self.f(x)
x[il] = xilold - x_eps
x[ir] = xirold + x_eps
f1m1 = self.f(x)
# x[il] = xilold - x_eps
x[ir] = xirold - x_eps
fm1m1 = self.f(x)
x[il] = xilold - 3. * x_eps
x[ir] = xirold + x_eps
f1m3 = self.f(x)
# x[il] = xilold - 3. * x_eps
x[ir] = xirold - x_eps
fm1m3 = self.f(x)
y = 0.0625 * (f13 - fm13 - 3. * (f11 -
fm11) + 3. * (f1m1 - fm1m1) - (f1m3 -
fm1m3)) / x_eps**4
elif complexity == 4:
xiold = x[i]
xjold = x[j]
xkold = x[k]
xlold = x[l]
x[i] = xiold + x_eps
x[j] = xjold + x_eps
x[k] = xkold + x_eps
x[l] = xlold + x_eps
f1111 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold + x_eps
# x[k] = xkold + x_eps
x[l] = xlold - x_eps
fm1111 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold + x_eps
x[k] = xkold - x_eps
x[l] = xlold + x_eps
f1m111 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold + x_eps
# x[k] = xkold - x_eps
x[l] = xlold - x_eps
fm1m111 = self.f(x)
# x[i] = xiold + x_eps
x[j] = xjold - x_eps
x[k] = xkold + x_eps
x[l] = xlold + x_eps
f11m11 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold - x_eps
# x[k] = xkold + x_eps
x[l] = xlold - x_eps
fm11m11 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold - x_eps
x[k] = xkold - x_eps
x[l] = xlold + x_eps
f1m1m11 = self.f(x)
# x[i] = xiold + x_eps
# x[j] = xjold - x_eps
# x[k] = xkold - x_eps
x[l] = xlold - x_eps
fm1m1m11 = self.f(x)
x[i] = xiold - x_eps
x[j] = xjold + x_eps
x[k] = xkold + x_eps
x[l] = xlold + x_eps
f111m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold + x_eps
# x[k] = xkold + x_eps
x[l] = xlold - x_eps
fm111m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold + x_eps
x[k] = xkold - x_eps
x[l] = xlold + x_eps
f1m11m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold + x_eps
# x[k] = xkold - x_eps
x[l] = xlold - x_eps
fm1m11m1 = self.f(x)
# x[i] = xiold - x_eps
x[j] = xjold - x_eps
x[k] = xkold + x_eps
x[l] = xlold + x_eps
f11m1m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold - x_eps
# x[k] = xkold + x_eps
x[l] = xlold - x_eps
fm11m1m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold - x_eps
x[k] = xkold - x_eps
x[l] = xlold + x_eps
f1m1m1m1 = self.f(x)
# x[i] = xiold - x_eps
# x[j] = xjold - x_eps
# x[k] = xkold - x_eps
x[l] = xlold - x_eps
fm1m1m1m1 = self.f(x)
y = 0.0625 * ((f1111 - fm1111) - (f1m111 -
fm1m111) - (f11m11 - fm11m11) + (f1m1m11 -
fm1m1m11) - (f111m1 - fm111m1) + (f1m11m1 -
fm1m11m1) + (f11m1m1 - fm11m1m1) - (f1m1m1m1 -
fm1m1m1m1)) / x_eps**4
elif complexity == 5:
xilold = x[il]
ximold = x[im]
xirold = x[ir]
x[il] = xilold + 2. * x_eps
x[im] = ximold + x_eps
x[ir] = xirold + x_eps
f112 = self.f(x)
# x[il] = xilold + 2 * x_eps
# x[im] = ximold + x_eps
x[ir] = xirold - x_eps
fm112 = self.f(x)
# x[il] = xilold + 2 * x_eps
x[im] = ximold - x_eps
x[ir] = xirold + x_eps
f1m12 = self.f(x)
# x[il] = xilold + 2 * x_eps
# x[im] = ximold - x_eps
x[ir] = xirold - x_eps
fm1m12 = self.f(x)
x[il] = xilold
x[im] = ximold + x_eps
x[ir] = xirold + x_eps
f110 = self.f(x)
# x[il] = xilold
# x[im] = ximold + x_eps
x[ir] = xirold - x_eps
fm110 = self.f(x)
# x[il] = xilold
x[im] = ximold - x_eps
x[ir] = xirold + x_eps
f1m10 = self.f(x)
# x[il] = xilold
# x[im] = ximold - x_eps
x[ir] = xirold - x_eps
fm1m10 = self.f(x)
x[il] = xilold - 2. * x_eps
x[im] = ximold + x_eps
x[ir] = xirold + x_eps
f11m2 = self.f(x)
# x[il] = xilold - 2. * x_eps
# x[im] = ximold + x_eps
x[ir] = xirold - x_eps
fm11m2 = self.f(x)
# x[il] = xilold - 2. * x_eps
x[im] = ximold - x_eps
x[ir] = xirold + x_eps
f1m1m2 = self.f(x)
# x[il] = xilold - 2. * x_eps
# x[im] = ximold - x_eps
x[ir] = xirold - x_eps
fm1m1m2 = self.f(x)
y = 0.0625 * (f112 - fm112 - (f1m12 - fm1m12)
- 2. * (f110 - fm110) + 2. * (f1m10 - fm1m10)
+ (f11m2 - fm11m2) - (f1m1m2 - fm1m1m2)) / x_eps**4
return y