"""Nonlinear equations.
We provide a variety of non-linear equations used for testing
numerical optimization algorithms.
"""
import functools as ft
import sys
import numdifftools as nd
import numpy as np
import pandas as pd
def _check_if_number(a, name):
r"""Function to check if object `a` is a number.
Parameters
----------
a : object
Object for which it should be tested whether it is a number or not.
name : str
String including the name of the object that appears in the error notification.
Examples
--------
>>> import numpy as np
>>> import pandas as pd
>>> import sys
>>>
>>> from temfpy.nonlinear_equations import _check_if_number
>>> a = 10
>>> _check_if_number(a, 'a')
True
"""
if not (isinstance(a, (int, float)) and not isinstance(a, bool)):
sys.exit(f"The parameter `{name}` must either be of type int or float.")
return True
def _check_if_val_x(x, name, length=None, length_type="equal"):
r"""Function to check if object `x` is an array, integer or float
and optionally if it has a specified length.
Currently considered as arrays are tuples, list, numpy.array
or pandas.Series type of objects.
Parameters
----------
x : object
Object for which it should be tested whether it is an array,
integer or float.
name : str
String including the name of the object that
appears in the error notification.
length : int
Desired length of array `x` for which the function
does not yield an error notification.
length_type : str
Indicates how the array length should be evaluated.
'equal' yields an error if the length of `x`
is not equal to the integer specified for `length`
and 'grtr_equ' yields an error if the length of `x`
is smaller.
Examples
--------
>>> import numpy as np
>>> import pandas as pd
>>> import sys
>>> from temfpy.nonlinear_equations import _check_if_val_x
>>>
>>> x = [10, 1, 4, 6]
>>> _check_if_val_x(x, 'x')
"""
if not isinstance(x, (int, float, list, tuple, pd.core.series.Series, np.ndarray)):
sys.exit(
f"The parameter `{name}` must either be an integer, float, "
f"list, numpy.array or pandas.Series.",
)
if not (isinstance(x, (int, float)) and not isinstance(x, bool)):
if length is not None:
if (length_type == "equal") and (len(x) != length):
sys.exit(f"The array `{name}` must have length {length}.")
if (length_type == "grtr_equ") and (len(x) < length):
sys.exit(f"The array `{name}` must have at least length {length}.")
def _exponential_val(x, a=10, b=1):
r"""Exponential function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_1(x) &= e^{x_1} - b \\
F_i(x) &= \frac{i}{a} (e^{x_i} +x_{i-1}) - b, i = 2,3, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
a : float, optional
The default value is 10.
b : float, optional
The default value is 1.
Returns
-------
array_like
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _exponential_val
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.normal(size = p)
>>> value = _exponential_val(x)
"""
_check_if_number(a, "a")
_check_if_number(b, "b")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
if isinstance(x, (int, float)):
p = 1
else:
p = len(x)
x = np.array(x)
x_im1 = np.concatenate((0, np.delete(x, p - 1)), axis=None)
rslt = (
(np.exp(x) + x_im1 - b)
* np.concatenate((a, np.array(range(2, p + 1))), axis=None)
/ a
)
return np.array(rslt)
def _exponential_jacobian(x, a=10):
r"""Analytical and numerical jacobian of the exponential function.
.. math::
F_1(x) &= e^{x_1} - b \\
F_i(x) &= \frac{i}{a} (e^{x_i} +x_{i-1}) - b, i = 2,3, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
a : float, optional
The default value is 10.
Returns
-------
numpy.array
Analytically derived Jacobian
numpy.array
Numerically derived Jacobian.
Only if dimension :math:`p > 1`.
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd #noqa
>>> from temfpy.nonlinear_equations import _exponential_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.normal(size = p)
>>> analytical_jacobian, numerical_jacobian = _exponential_jacobian(x)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_number(a, "a")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
if isinstance(x, (int, float)):
p = 1
else:
p = len(x)
x = np.array(x)
diag_mat = np.diag(
np.exp(np.array(x))
* np.array(range(1, p + 1))
/ np.append([1], np.repeat(a, p - 1)),
)
off_diag_mat = np.diag(
np.array(range(2, p + 1)) / np.array(np.repeat(a, p - 1)),
k=-1,
)
jacobian = diag_mat + off_diag_mat
j_numdiff = nd.Jacobian(_exponential_val)
jacobian_numdiff = j_numdiff(np.array(x))
return jacobian, jacobian_numdiff
[docs]def exponential(x, a=10, b=1, jac=False):
r"""Exponential function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_1(x) &= e^{x_1} - b \\
F_i(x) &= \frac{i}{a} (e^{x_i} +x_{i-1}) - b \\
& \quad i = 2,3, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
a : float, optional
The default value is 10.
b : float, optional
The default value is 1.
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian and the
numerically derived Jacobian.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import exponential
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(1,20)
>>> x = np.zeros(p)
>>> np.allclose(exponential(x), np.zeros(p))
True
"""
if a == 0:
sys.exit("a must be different from 0.")
if jac is False:
return _exponential_val(x, a=a, b=b)
if jac is True:
return _exponential_val(x, a=a, b=b), _exponential_jacobian(x, a=a)
def _trig_exp_i(xi, a=(3, 2, 5, 4, 3, 2, 8, 4, 3)):
r"""Trigonometrical exponential function. Used to build
the function :func:`trig_exp_val`.
.. math::
F_i(x) = - x_{i-1}e^(x_{i-1} - x_i) + x_i(a_4+a_5x_i^2)
+ a_6x_{i+1} + \sin(x_i - x_{i+1})\sin(x_i + x_{i+1}) - a_7,
i = 2,3, \dots, p-1
Parameters
----------
x_i : array_like
Input domain with dimension :math:`3`.
a : array_like, optional
The default array is (3,2,5,4,3,2,8,4,3).
Returns
-------
float
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _trig_exp_i
>>>
>>> p = 3
>>> np.random.seed(123)
>>> x = np.random.normal(size = p)
>>> F_i = _trig_exp_i(x)
"""
_check_if_val_x(a, "a", length=9)
_check_if_val_x(xi, "xi", length=3)
xi = np.array(xi)
a = np.array(a)
rslt = (
-xi[0] * np.exp(xi[0] - xi[1])
+ xi[1] * (a[3] + a[4] * xi[1] ** 2)
+ a[5] * xi[2]
+ np.sin(xi[1] - xi[2]) * np.sin(xi[1] + xi[2])
- a[6]
)
return rslt
def _trig_exp_val(x, a=(3, 2, 5, 4, 3, 2, 8, 4, 3)):
r"""Trigonometrical exponential function.
.. math::
F_1(x) &= a_1x_1^3 + a_2x_2 - a_3 + \sin(x_1 - x_2)\sin(x1+x2) \\
F_i(x) &= - x_{i-1}e^{x_{i-1} - x_i} + x_i(a_4+a_5x_i^2)
+ a_6x_{i+1} + \sin(x_i - x_{i+1})\sin(x_i + x_{i+1}) - a_7,
i = 2,3, \dots, p-1 \\
F_p(x) &= -x_{p-1}e^{x_{p-1}-x_p} + a_8x_p - a_9
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
a : array_like, optional
The default array is [3,2,5,4,3,2,8,4,3].
Returns
-------
array_like
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _trig_exp_val
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.normal(size = p)
>>> value = _trig_exp_val(x)
"""
_check_if_val_x(a, "a", length=9)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
rslt = [
a[0] * x[0] ** 3
+ a[1] * x[1]
- a[2]
+ np.sin(x[0] - x[1]) * np.sin(x[0] + x[1]),
]
for i in range(2, p):
rslt.append(_trig_exp_i(xi=x[(i - 2) : (i + 1)], a=a))
rslt.append(-x[p - 2] * np.exp(x[p - 2] - x[p - 1]) + a[7] * x[p - 1] - a[8])
return np.array(rslt)
def _trig_exp_jacobian(x, a=(3, 2, 5, 4, 3, 2, 8, 4, 3)):
r"""Trigonometrical exponential function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p \end{pmatrix}^T \\
F_1(x) &= a_1x_1^3 + a_2x_2 - a_3 + \sin(x_1 - x_2)\sin(x1+x2) \\
F_i(x) &= - x_{i-1}e^{x_{i-1} - x_i} + x_i(a_4+a_5x_i^2)
+ a_6x_{i+1} \\
& \quad + \sin(x_i - x_{i+1})\sin(x_i + x_{i+1}) - a_7,
i = 2,3, \dots, p-1 \\
F_p(x) &= -x_{p-1}e^{x_{p-1}-x_p} + a_8x_p - a_9
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
a : array_like, optional
The default array is (3,2,5,4,3,2,8,4,3).
Returns
-------
array_like
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _trig_exp_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.normal(size = p)
>>> analytical_jacobian, numerical_jacobian = _trig_exp_jacobian(x)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_val_x(a, "a", length=9)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
x_im1 = np.delete(x, [p - 1, p - 2])
x_i = np.delete(x, [0, p - 1])
x_ip1 = np.delete(x, [0, 1])
diag_mat = np.diag(
np.concatenate(
(
np.array(
3 * a[0] * x[0] ** 2
+ np.sin(x[0] - x[1]) * np.cos(x[0] + x[1])
+ np.cos(x[0] - x[1]) * np.sin(x[0] + x[1]),
),
(
x_im1 * np.exp(x_im1 - x_i)
+ np.repeat(a[3], p - 2)
+ 3 * np.repeat(a[4], p - 2) * x_i ** 2
+ np.sin(x_i - x_ip1) * np.cos(x_i + x_ip1)
+ np.cos(x_i - x_ip1) * np.sin(x_i + x_ip1)
),
np.array(x[p - 2] * np.exp(x[p - 2] - x[p - 1]) + a[7]),
),
axis=None,
),
)
off_diag_p1_mat = np.diag(
(
np.concatenate((a[1], np.repeat(a[5], p - 2)), axis=None)
+ np.sin(np.delete(x, [p - 1]) - np.delete(x, [0]))
* np.cos(np.delete(x, [p - 1]) + np.delete(x, [0]))
- np.cos(np.delete(x, [p - 1]) - np.delete(x, [0]))
* np.sin(np.delete(x, [p - 1]) + np.delete(x, [0]))
),
k=1,
)
off_diag_m1_mat = np.diag(
-(np.delete(x, [p - 1]) + np.repeat(1, p - 1))
* np.exp(np.delete(x, [p - 1]) - np.delete(x, [0])),
k=-1,
)
jacobian = off_diag_p1_mat + off_diag_m1_mat + diag_mat
j_numdiff = nd.Jacobian(_trig_exp_val)
jacobian_numdiff = j_numdiff(np.array(x))
return jacobian, jacobian_numdiff
[docs]def trig_exp(x, a=(3, 2, 5, 4, 3, 2, 8, 4, 3), jac=False):
r"""Trigonometrical exponential function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_1(x) &= a_1x_1^3 + a_2x_2 - a_3 + \sin(x_1 - x_2)\sin(x_1+x_2) \\
F_i(x) &= - x_{i-1}e^{x_{i-1} - x_i} + x_i(a_4+a_5x_i^2)
+ a_6x_{i+1} \\
& \quad + \sin(x_i - x_{i+1})\sin(x_i + x_{i+1}) - a_7 \\
& \quad i = 2,3, \dots, p-1 \\
F_p(x) &= -x_{p-1}e^{x_{p-1}-x_p} + a_8x_p - a_9
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
a : array_like, optional
The default array is (3,2,5,4,3,2,8,4,3).
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian
and the numerically derived Jacobian.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import trig_exp
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(3,20)
>>> x = np.zeros(p)
>>> compare = np.insert(np.array([-5,-3]), 1, np.repeat(-8, p-2))
>>> np.allclose(trig_exp(x), compare)
True
"""
if jac is False:
return _trig_exp_val(x, a=a)
else:
return _trig_exp_val(x, a=a), _trig_exp_jacobian(x, a=a)
def _broyden_val(x, a=(3, 0.5, 2, 1)):
r"""Broyden tridiagonal function.
.. math::
F_1(x) &= x_1(a_1 - a_2 x_1) -a_3 x_{2} + a_4 \\
F_i(x) &= x_i(a_1 - a_2 x_i)-x_{i-1} -a_3 x_{i+1}
+ a_4, i = 2,3, \dots, p-1 \\
F_p(x) &= x_p(a_1 - a_2 x_p)-x_{p-1} + a_4
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
a : array_like, optional
The default array is (3, 0.5, 2, 1)
Returns
-------
float
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _broyden_val
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = - np.random.uniform(size = p)
>>> value = _broyden_val(x)
"""
_check_if_val_x(a, "a", length=4)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
x_ip1 = np.concatenate((np.delete(x, 0), 0), axis=None)
x_im1 = np.concatenate((0, np.delete(x, p - 1)), axis=None)
rslt = x * (a[0] - a[1] * x) - x_im1 - a[2] * x_ip1 + 1
return np.array(rslt)
def _broyden_jacobian(x, a=(3, 0.5, 2, 1)):
r"""Broyden tridiagonal function.
.. math::
F_1(x) &= x_1(a_1 - a_2 x_1) -a_3 x_{2} + a_4 \\
F_i(x) &= x_i(a_1 - a_2 x_i)-x_{i-1} -a_3 x_{i+1}
+ a_4, i = 2,3, \dots, p-1 \\
F_p(x) &= x_p(a_1 - a_2 x_p)-x_{p-1} + a_4
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
a : array_like, optional
The default array is (3, 0.5, 2, 1)
Returns
-------
numpy.array
Analytically derived Jacobian
numpy.array
Numerically derived Jacobian
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _broyden_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = - np.random.uniform(size = p)
>>> analytical_jacobian, numerical_jacobian = _broyden_jacobian(x)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_val_x(a, "a", length=4)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
jacobian = (
np.diag(np.repeat(a[0], p) - x * np.repeat(2 * a[1], p))
+ np.diag(np.repeat(-a[2], p - 1), k=1)
+ np.diag(np.repeat(-1, p - 1), k=-1)
)
j_numdiff = nd.Jacobian(_broyden_val)
jacobian_numdiff = j_numdiff(np.array(x))
return jacobian, jacobian_numdiff
[docs]def broyden(x, a=(3, 0.5, 2, 1), jac=False):
r"""Broyden tridiagonal function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_1(x) &= x_1(a_1 - a_2 x_1) -a_3 x_{2} + a_4 \\
F_i(x) &= x_i(a_1 - a_2 x_i)-x_{i-1} -a_3 x_{i+1} + a_4 \\
& \quad i = 2,3, \dots, p-1 \\
F_p(x) &= x_p(a_1 - a_2 x_p)-x_{p-1} + a_4
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
a : array_like, optional
The default array is (3, 0.5, 2, 1).
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian
and the numerically derived Jacobian.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import broyden
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(3,20)
>>> x = np.zeros(p)
>>> np.allclose(broyden(x), np.repeat(1,p))
True
"""
if jac is False:
return _broyden_val(x=x, a=a)
else:
return _broyden_val(x=x, a=a), _broyden_jacobian(x=x, a=a)
def _rosenbrock_ext_val(x, a=(10, 1)):
r"""Extended-Rosenbrock function.
.. math::
F_{2i-1}(x) &= a_1(x_{2i} - x_{2i-1}^2) \\
F_{2i}(x) &= a_2 - x_{2i-1}, i = 1,2,3, \dots, \frac{p}{2}
Parameters
----------
x : array_like
Input domain with even dimension :math:`p > 1`.
a : array_like, optional
The default array is (10,1)
Returns
-------
float
Output domain
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _rosenbrock_ext_val
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = - np.random.uniform(size = p)
>>> value = _rosenbrock_ext_val(x)
"""
_check_if_val_x(a, "a", length=2)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
if p % 2 != 0:
sys.exit("x must consist of an even number of parameters.")
xl = np.concatenate((np.delete(x, 0), 0), axis=None)
xh = np.concatenate((np.delete(x, p - 1), 0), axis=None)
f_odd = a[0] * (xl - xh ** 2) * np.resize((1, 0), p)
f_even = np.delete(1 - np.concatenate((0, x), axis=None), p) * np.resize((0, 1), p)
rslt = f_odd + f_even
return np.array(rslt)
def _rosenbrock_ext_jacobian(x, a=(10, 1)):
r"""Extended-Rosenbrock function.
.. math::
F_{2i-1}(x) &= a_1(x_{2i} - x_{2i-1}^2) \\
F_{2i}(x) &= a_2 - x_{2i-1}, i = 1,2,3, \dots, \frac{p}{2}
Parameters
----------
x : array_like
Input domain with even dimension :math:`p > 1`.
a : array_like, optional
The default array is (10,1)
Returns
-------
numpy.array
Analytically derived Jacobian
numpy.array
Numerically derived Jacobian
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _rosenbrock_ext_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.uniform(size = p)
>>> analytical_jacobian, numerical_jacobian = _rosenbrock_ext_jacobian(x)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_val_x(a, "a", length=2)
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
a = np.array(a)
p = len(x)
if p % 2 != 0:
sys.exit("x must consist of an even number of parameters.")
diag_mat = np.diag(np.repeat(-2, p) * np.repeat(a[0], p) * x * np.resize([1, 0], p))
off_diag_p1_mat = np.diag(np.resize([a[0], 0], p - 1), k=1)
off_diag_m1_mat = np.diag(np.resize([-1, 0], p - 1), k=-1)
jacobian = diag_mat + off_diag_p1_mat + off_diag_m1_mat
j_numdiff = nd.Jacobian(_rosenbrock_ext_val)
jacobian_numdiff = j_numdiff(np.array(x))
return jacobian, jacobian_numdiff
[docs]def rosenbrock_ext(x, a=(10, 1), jac=False):
r"""Extended-Rosenbrock function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_{2i-1}(x) &= a_1(x_{2i} - x_{2i-1}^2) \\
F_{2i}(x) &= a_2 - x_{2i-1}, \\
& \quad i = 1,2,3, \dots, \frac{p}{2}
Parameters
----------
x : array_like
Input domain with even dimension :math:`p > 1`.
a : array_like, optional
The default array is (10,1).
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian and the
numerically derived Jacobian.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import rosenbrock_ext
>>>
>>> np.random.seed(123)
>>> p = 2*np.random.randint(1,20)
>>> x = np.zeros(p)
>>> compare = np.resize([0,1], p)
>>> np.allclose(rosenbrock_ext(x), compare)
True
"""
if jac is False:
return _rosenbrock_ext_val(x=x, a=a)
else:
return _rosenbrock_ext_val(x=x, a=a), _rosenbrock_ext_jacobian(x=x, a=a)
def _troesch_val(x, rho=10, a=2):
r"""Troesch function.
.. math::
F_1(x) &= a_1x_1 + \rho h^2 \sinh(\rho x_1) - x_{2}, \\
F_i(x) &= a_1x_i + \rho h^2 \sinh(\rho x_i) - x_{i-1} - x_{i+1}, i = 2,3,
\dots, p-1\\
F_p(x) &= a_1x_p + \rho h^2 \sinh(\rho x_p) - x_{p-1}
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
rho : float, optional
The default value is 10
a : float, optional
The default value is 2
Returns
-------
array_like
Output domain
Notes
-----
:math:'h = \frac{1}{p+1}'
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _troesch_val
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.uniform(size = p)
>>> val = _troesch_val(x)
"""
_check_if_number(rho, "rho")
_check_if_number(a, "a")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
p = len(x)
h = 1 / (p + 1)
x_ip1 = np.concatenate((np.delete(x, 0), 0), axis=None)
x_im1 = np.concatenate((0, np.delete(x, p - 1)), axis=None)
rslt = a * x + rho * h ** 2 * np.sinh(rho * x) - x_ip1 - x_im1
return np.array(rslt)
def _troesch_jacobian(x, rho=10, a=2):
r"""Troesch function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
h &= \frac{1}{p+1} \\
F_1(x) &= ax_1 + \rho h^2 \sinh(\rho x_1) - x_{2}, \\
F_i(x) &= ax_i + \rho h^2 \sinh(\rho x_i) - x_{i-1} - x_{i+1} \\
& \quad i = 2,3, \dots, p-1 \\
F_p(x) &= ax_p + \rho h^2 \sinh(\rho x_p) - x_{p-1
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
rho : float, optional
The default value is 10
a : float, optional
The default value is 2
Returns
-------
numpy.array
Analytically derived Jacobian
numpy.array
Numerically derived Jacobian
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _troesch_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.uniform(size = p)
>>> analytical_jacobian, numerical_jacobian = _troesch_jacobian(x)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_number(rho, "rho")
_check_if_number(a, "a")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
x = np.array(x)
p = len(x)
h = 1 / (p + 1)
diag_mat = np.diag(
np.repeat(a, p)
+ np.repeat(rho ** 2, p)
* np.repeat(h ** 2, p)
* np.cosh(np.repeat(rho, p) * x),
)
off_diag_p1_mat = np.diag(np.repeat(-1, p - 1), k=1)
off_diag_m1_mat = np.diag(np.repeat(-1, p - 1), k=-1)
jacobian = diag_mat + off_diag_p1_mat + off_diag_m1_mat
j_numdiff = nd.Jacobian(_troesch_val)
jacobian_numdiff = j_numdiff(np.array(x))
return jacobian, jacobian_numdiff
[docs]def troesch(x, rho=10, a=2, jac=False):
r"""Troesch function.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
h &= \frac{1}{p+1} \\
F_1(x) &= ax_1 + \rho h^2 \sinh(\rho x_1) - x_{2}, \\
F_i(x) &= ax_i + \rho h^2 \sinh(\rho x_i) - x_{i-1} - x_{i+1} \\
& \quad i = 2,3, \dots, p-1 \\
F_p(x) &= ax_p + \rho h^2 \sinh(\rho x_p) - x_{p-1}
Parameters
----------
x : array_like
Input domain with dimension :math:`p > 1`.
rho : float, optional
The default value is 10.
a : float, optional
The default value is 2.
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian and the
numerically derived Jacobian.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import troesch
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(1,20)
>>> x = np.zeros(p)
>>> np.allclose(troesch(x), np.zeros(p))
True
"""
if jac is False:
return _troesch_val(x=x, rho=rho, a=a)
else:
return _troesch_val(x=x, rho=rho, a=a), _troesch_jacobian(x=x, rho=rho, a=a)
def _chandrasekhar_val(x, y, c):
r"""Discretized version of Chandrasekhar’s H-equation.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_i(x) &= x_i - \left(1 - \frac{c}{2p} \sum^p_{j=1}
\frac{y_i x_j}{y_i + y_j} \right)^{-1} \\
& \quad i = 1,2, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
y : array_like,
Array of constants with dimension :math:`p`
c : float
Constant parameter
Returns
-------
array_like
Output domain
array_like
Tuple containing the analytically derived Jacobian and the
numerically derived Jacobian
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _chandrasekhar_val
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(1,20)
>>> x = np.repeat(2,p)
>>> y = np.repeat(1,p)
>>> c = 1
>>> np.allclose(_chandrasekhar_val(x,y,c), np.zeros(p))
True
"""
_check_if_number(c, "c")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
_check_if_val_x(y, "y")
if isinstance(x, (int, float)) and isinstance(y, (int, float)):
p = 1
term_sum = x / (2 * y)
else:
p = len(x)
if len(x) != len(y):
sys.exit("The arrays `x` and `y` must have the same length.")
x = np.array(x)
y = np.array(y)
term_sum = []
for i in range(0, p):
term_sum.append(np.sum(x / (y[i] + y)))
rslt = x - 1 / (1 - c * y / (2 * p) * term_sum)
return np.array(rslt)
def _chandrasekhar_jacobian(x, y, c):
r"""Discretized version of Chandrasekhar’s H-equation.
.. math::
F_i(x) = x_i - \left(1 - \frac{c}{2p} \sum^p_{j=1}
\frac{y_i x_j}{y_i + y_j} \right)^{-1}
i = 1,2, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
y : array_like,
Array of constants with dimension :math:'p'
c : float
Constant parameter
Returns
-------
numpy.array
Analytically derived Jacobian
numpy.array
Numerically derived Jacobian.
Only if dimension :math:`p > 1`.
Examples
--------
>>> import numpy as np
>>> import numdifftools as nd
>>> from temfpy.nonlinear_equations import _chandrasekhar_jacobian
>>>
>>> p = 10
>>> np.random.seed(123)
>>> x = np.random.uniform(size = p)
>>> y = np.random.normal(size = p)
>>> c = 2
>>> analytical_jacobian, numerical_jacobian = _chandrasekhar_jacobian(x,y, c)
>>> np.allclose(analytical_jacobian, numerical_jacobian)
True
"""
_check_if_number(c, "c")
_check_if_val_x(x, "x", length=1, length_type="grtr_equ")
_check_if_val_x(y, "y")
if isinstance(x, (int, float)) and isinstance(y, (int, float)):
p = 1
jacobian = (
-1 / (1 - c * y / (2 * p) * x / (2 * y)) ** 2 * (c * y / (2 * p) / (2 * y))
)
else:
p = len(x)
if len(x) != len(y):
sys.exit("The arrays `x` and `y` must have the same length.")
x = np.array(x)
y = np.array(y)
jacobian = np.zeros((p, p))
for k in range(0, p):
term_sum = []
for i in range(0, p):
term_sum.append(np.sum(x / (y[i] + y)))
column_k = (
-1
/ (1 - c * y / (2 * p) * term_sum) ** 2
* (c * y / (2 * p) / (y + y[k]))
)
jacobian[:, k] = column_k
jacobian = jacobian + np.identity(p)
_chandrasekhar_help_num_ft = ft.partial(_chandrasekhar_val, y=y, c=c)
j_numdiff = nd.Jacobian(_chandrasekhar_help_num_ft)
jacobian_numdiff = j_numdiff(x)
# return more dimensional cas
return jacobian, jacobian_numdiff
# return one dimensional case
return jacobian
[docs]def chandrasekhar(x, y, c, jac=False):
r"""Discretized version of Chandrasekhar’s H-equation.
.. math::
x &\mapsto \begin{pmatrix} F_1(x) & F_2(x) & \dots & F_p(x) \end{pmatrix}^T \\
F_i(x) &= x_i - \left(1 - \frac{c}{2p} \sum^p_{j=1}
\frac{y_i x_j}{y_i + y_j} \right)^{-1} \\
& \quad i = 1,2, \dots, p
Parameters
----------
x : array_like
Input domain with dimension :math:`p`.
y : array_like,
Array of constants with dimension :math:`p`.
c : float
Constant parameter.
jac : bool
If True, an additional array containing the numerically derived Jacobian
and the analytically derived Jacobian is returned.
The default is False.
Returns
-------
array_like
Output domain
array_like
Only returned if :math:`jac = True`.
Tuple containing the analytically derived Jacobian and the
numerically derived Jacobian.
Numerically derived Jacobian returned only if dimension :math:`p > 1`.
References
----------
.. [V2009] Varadhan, R., and Gilbert, P. D. (2009). BB: An R package for solving a
large system of nonlinear equations and for optimizing a high-dimensional
nonlinear objective function. *Journal of Statistical Software*,
32(1): 1–26.
Examples
--------
>>> import numpy as np
>>> from temfpy.nonlinear_equations import chandrasekhar
>>>
>>> np.random.seed(123)
>>> p = np.random.randint(1,20)
>>> x = np.repeat(2,p)
>>> y = np.repeat(1,p)
>>> c = 1
>>> np.allclose(chandrasekhar(x,y,c), np.zeros(p))
True
"""
x = np.atleast_1d(x)
p = len(x)
x_check_zero = np.matrix(x).T
x_matrix = x_check_zero @ np.matrix(np.repeat(1, p))
matrix_shouldnt_have_zeros = x_matrix + x_matrix.T
if 0 in matrix_shouldnt_have_zeros:
sys.exit("Every sum of two numbers from y must be different from 0.")
if jac is False:
return _chandrasekhar_val(x=x, y=y, c=c)
else:
return _chandrasekhar_val(x=x, y=y, c=c), _chandrasekhar_jacobian(x=x, y=y, c=c)