"""
nle module description
"""
import numpy as np
from scipy import optimize
from nutopy.mod_nlesolve import mod_nlesolve
from nutopy.errors import *
from nutopy import options
from nutopy.tools import *
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
def finite_diff(fun, x, *args, **kwargs):
v_eps = np.finfo(float).eps
if isinstance(x, float):
t = np.sqrt(v_eps) * np.sqrt(np.maximum(1.0, np.abs(x)))
j = (fun(x + t, *args, **kwargs) - fun(x, *args, **kwargs)) / t
else:
n = x.size
j = np.zeros((n, n))
f = fun(x, *args, **kwargs)
for i in range(0, n):
t = np.sqrt(v_eps) * np.sqrt(np.maximum(1.0, np.abs(x[i])))
xi = x[i]
x[i] = xi + t
j[:, i] = (fun(x, *args, **kwargs) - f) / t
x[i] = xi
return j
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def solve(f, x0, *, df=None, options=None, callback=None, args=(), kwargs=None):
"""
Solver of nonlinear equations:
.. math::
f(x) = 0,
with :math:`f : \mathbf{R}^n \\rightarrow \mathbf{R}^n`.
The solver is iterative and so it constructs a sequence of iterates :math:`x^{(k)}`
which must converge to a solution denoted :math:`x^*` of :math:`f(x)=0`. The user
must provide an ``initial guess`` :math:`x^{(0)}`.
Parameters
----------
f : callable
A vector function representing ``f(x)``.
x0 : float vector or :class:`float`
Initial guess :math:`x^{(0)}`.
df : callable, optional, key only
The Jacobian ``f'(x)`` which returns a float matrix.
options : :meth:`nutopy.nle.Options`, optional, key only
A dictionnary of solver options.
callback: callable, optional, key only
A function of the form ``callback(infos)`` which returns a flag
(named ``status``).
The ``nle.solve`` method stops if ``SolverMethod=='hybrj'`` and ``status < 0``
at some iteration.
The argument ``infos`` is a structure containing:
- ``infos.x`` : ``x`` at current iteration;
- ``infos.f`` : ``f(x)`` at current iteration.
args : tuple, key only
optional arguments for ``f`` (and ``df`` if provided).
kwargs : dictionnary, key only
keywords optional arguments for ``f`` (and ``df`` if provided).
Returns
-------
sol
A dictionary/struct of outputs:
- ``sol.x`` or ``sol.get('x')`` (float vector) : solution of f(x) = 0 if
the method has converged;
- ``sol.f`` (float vector) : value ``f(xsol)``;
- ``sol.nfev`` (:class:`int`) : number of evaluations of ``f``;
- ``sol.njev`` (:class:`int`) : number of evaluations of ``df``;
- ``sol.status`` (:class:`int`) : termination status of the solver;
- ``sol.success`` (:class:`bool`) : whether or not the solver exited successfully;
- ``sol.message`` (:class:`str`) : description of the cause of the termination.
See Also
--------
nutopy.nle.Options
nutopy.errors
Notes
-----
- Use :class:`numpy.ndarray` to define float vectors and matrices.
- You can define a unique function ``f`` coding f and its derivative:
.. highlight:: python
>>> from nutopy.tools import *
>>> def df(x, dx):
... dy = 2.0 * dx
... return dy
>>> @tensorize(df)
... def f(x):
... y = 2.0 * x - 1.0
... return y
Then, call ``nle.solve`` with ``df=f``:
.. highlight:: python
>>> from nutopy import nle
>>> nle.solve(f, x0=1.0, df=f)
>>> # or equivalently: df is not required if f is tensorized
>>> nle.solve(f, x0=1.0)
Examples
--------
>>> import numpy as np
>>> from nutopy import nle
Example 1: scalar
>>> sol = nle.solve(lambda x: 2.0 * x - 1.0, 0.0)
>>> sol.x
0.5
Example 2: use numpy array
>>> def f(x):
... y = np.zeros(x.size)
... y[0] = 2.0*x[0]-1.0
... y[1] = np.exp(x[1])-1.0
... return y
>>> x0 = np.array([10.0, 2.0])
>>> sol = nle.solve(f, x0)
Example 3: with Jacobian and options
>>> def df(x):
... j = np.eye(x.size)
... j[0,0] = 2.0
... j[1,1] = np.exp(x[1])
... return j
>>> opt = nle.Options(SolverMethod='hybrj', Display='on', TolX=1e-8)
>>> sol = nle.solve(f, x0, df=df, options=opt)
Example 4: with optional arguments
>>> def f(x, a, b, c=1.0):
... y = np.zeros(x.size)
... y[0] = a*x[0]-b
... y[1] = np.exp(x[1])-c
... return y
>>> a = 2.0
>>> b = 1.0
>>> x0 = np.array([10.0, 2.0])
>>> sol = nle.solve(f, x0, args=(a, b), kwargs={'c' : 0.1})
# Example 5: with callback
>>> def f(x):
... y = np.zeros(x.size)
... y[0] = 2.0 * x[0] - 1.0
... return y
>>> def cb(infos):
... return 0
>>> x0 = np.array([0.0])
>>> sol = nle.solve(f, x0, callback=cb)
# Example 6: with df=f
>>> from nutopy.tools import *
>>> def df(x):
... dy = np.array([2.0 * dx[0], np.exp(x[1]) * dx[1]])
... return dy
>>> @tensorize(df)
... def df(x):
... y = np.array([2.0 * x[0] - 1.0, np.exp(x[1]) - 1.0])
... return y
>>> sol = nle.solve(f, np.array([10.0, 2.0]), df=f)
Raises
------
ArgumentDimensionError
Variable x0 must be a one dimensional numpy array or a float.
ArgumentTypeError
Variable options must be a Options object.
"""
# Cast of input variables
if kwargs is None:
kwargs = {}
if not isinstance(x0, float):
x0 = np.asarray(x0)
# Tests on input variables
if options is None:
options = Options()
# Variable x0 must be a vector or a float
msg_erreur = 'Variable x0 must be a one dimensional numpy array or a float.';
if isinstance(x0, np.ndarray):
# Variable x0 must be a vector
if x0.ndim == 1:
n = x0.size
elif x0.ndim == 2:
(lig, col) = x0.shape
if lig == 1 or col == 1:
n = x0.size
else:
raise ArgumentDimensionError(msg_erreur)
else:
raise ArgumentDimensionError(msg_erreur)
elif isinstance(x0, float):
n = 1
else:
raise ArgumentTypeError(msg_erreur)
# Variable options must be a Options object
if not isinstance(options, Options):
raise ArgumentTypeError('Variable options must be a Options object.')
# create callback function passed to fortran routine
display = options.get('Display')=='on'
if not isinstance(args,tuple):
args = (args,)
#
if f is df:
jac_case = 2 # f is tensorized
elif df is None:
#if istensorized(f, (1,), x0, *args, **kwargs):
# df = f
# jac_case = 2 # f is tensorized
#else:
jac_case = 0 # no jacobian => finite diff
else:
jac_case = 1 # the Jacobian is given
#
if jac_case == 0: # No jacobian
myjac = lambda x: finite_diff(f, x, *args, **kwargs) # finite diff
elif jac_case == 2: # f is tensorized
def myjac(x):
if isinstance(x, float):
_, dy = df((x, 1.0), *args, **kwargs)
else:
n = x.size
dy = np.zeros((n,n))
Id = np.eye(n)
for i in range(0, n): _, dy[:,i] = df((x, Id[:, i]), *args, **kwargs)
return dy
else: # jac_case = 1: the Jacobian is given
myjac = lambda x: df(x, *args, **kwargs)
myfun = lambda x: f(x, *args, **kwargs)
if callback is None:
def callback_nle(x, f):
if display:
callbackPrint(x, f)
return 0
else:
def callback_nle(x, f):
if display:
callbackPrint(x, f)
infos_dict = dict(x=x, f=f)
infos = NLEInfos(infos_dict)
if callback.__code__.co_argcount > 1:
return callback(infos, *args, **kwargs)
else:
return callback(infos) # we can call callback without optional
# arguments
# prepare call to Fortran solver
[dw, iw, lsw, sw] = options.getAllOptions()
# Code ascii de la chaine de caractere
solver_method = options.get('SolverMethod')
if solver_method == 'hybrj':
if not isinstance(x0, float):
x0.astype(dtype='d', order='F')
[xsol, fsol, nfev, njev, flag, messageFortran] = mod_nlesolve.nlesolvepy( \
myfun, myjac, x0, np.array([1.0]), callback_nle, \
dw, iw, [ord(c) for c in sw], lsw, n, 1, len(dw), len(iw), len(lsw), len(sw))
else:
[xsol, fsol, nfev, njev, flag, messageFortran] = mod_nlesolve.nlesolvepy_float( \
myfun, myjac, x0, callback_nle, \
dw, iw, [ord(c) for c in sw], lsw, len(dw), len(iw), len(lsw), len(sw))
message = messageFortran.strip().decode()
status = flag
if status == 1:
success = True
else:
success = False
else:
if solver_method == 'hybr' or solver_method == 'lm':
sol = optimize.root(myfun, x0, jac=myjac, method=solver_method, callback=None,
options=options.get('ScipyRootOptions'))
else:
if callback is None:
sol = optimize.root(myfun, x0, jac=None, method=solverMethod, \
callback=callbackPrint, options=options.get('ScipyRootOptions'))
else:
sol = optimize.root(myfun, x0, jac=None, method=solver_method, \
callback=callback_nle, options=options.get('ScipyRootOptions'))
xsol = sol.x
fsol = sol.fun
message = sol.message
success = sol.success
status = sol.status
try:
nfev = sol.nfev
except AttributeError as err:
nfev = -1
try:
njev = sol.njev
except AttributeError as err:
njev = -1
# print('\n Results of the nle solver method:\n')
# print(' xsol = ', xsol, '\n');
# print(' f(xsol) = ', fsol, '\n');
# print(' nfev = ', nfev, '\n');
# print(' njev = ', njev, '\n');
# print(' status = ', status, '\n');
# print(message)
# post computations, prepare outputs
#if n == 1:
if isinstance(x0, float):
xsol = type(x0)(xsol)
fsol = type(x0)(fsol)
# -----------------------------------------------------------------------------------
sol = dict(x=xsol, f=fsol, nfev=nfev, njev=njev, status=status, \
success=success, message=message)
nlesol = NLESol(sol)
# -----------------------------------------------------------------------------------
# display outputs
if display:
print('\n Results of the nle solver method:\n')
print(' xsol = ', nlesol.x);
print(' f(xsol) = ', nlesol.f);
print(' nfev = ', nlesol.nfev);
print(' njev = ', nlesol.njev);
print(' status = ', nlesol.status);
print(' success = ', nlesol.success, '\n');
print(' ' + nlesol.message + '\n')
# -----------------------------------------------------------------------------------
# reset counter
callbackPrint.counter=0
return nlesol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
def static_vars(**kwargs):
def decorate(func):
for k in kwargs:
setattr(func, k, kwargs[k])
return func
return decorate
@static_vars(counter=0)
def callbackPrint(x,f):
if(callbackPrint.counter==0):
print('\n Calls |f(x)| |x|\n ')
print('{0:10}'.format(callbackPrint.counter+1) + \
'{0:23.15e}'.format(np.linalg.norm(f)) + \
'{0:23.15e}'.format(np.linalg.norm(x)))
callbackPrint.counter += 1
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class NLESol:
def __init__(self, *args, **kwargs):
"""
Parameters
----------
args : tuple
Anonymous arguments
kwargs :
Dictionary of named arguments
"""
self.sol = dict()
self.update(*args, **kwargs)
self.x = self.sol['x']
self.f = self.sol['f']
self.nfev = self.sol['nfev']
self.njev = self.sol['njev']
self.success = self.sol['success']
self.message = self.sol['message']
self.status = self.sol['status']
def update(self, *args, **kwargs):
"""
Doc de update
"""
for dico in args:
if isinstance(dico, dict):
self.__private_update(dico)
else:
raise ArgumentTypeError(
'NLESol requires arguments of type dict or a list of param=value, \
see help(NLESol)!')
self.__private_update(kwargs)
def __private_update(self, dictionary):
"""
Doc de __private_update
"""
for key, val in dictionary.items():
if key in self.sol:
if isinstance(self.sol[key], val.__class__):
self.sol[key] = val
else:
raise ArgumentTypeError('NLESol ' + key + ' must be of type: ' + \
str(self.sol[key].__class__) + \
' instead of type: ' + str(val.__class__) + '!')
else:
self.sol[key] = val
def get(self, key):
"""
Doc de get
"""
try:
value = self.sol[key]
except KeyError as err:
raise OptionNameError(key + ' is not a NLESol, see help(NLESol)!\n')
return value
def __str__(self):
"""
Doc de __str__
"""
keys = sorted(self.sol.keys())
string = self.__class__.__name__ + ':\n\n'
for key in keys:
string = string + "\t" + key.ljust(20) + ' = ' + str(self.sol[key]) + '\n'
return string
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class NLEInfos:
def __init__(self, *args, **kwargs):
"""
Parameters
----------
args : tuple
Anonymous arguments
kwargs : dictionary
Named arguments
"""
self.sol = dict()
self.update(*args, **kwargs)
self.x = self.sol['x']
self.f = self.sol['f']
def update(self, *args, **kwargs):
"""
Doc de update
"""
for dico in args:
if isinstance(dico, dict):
self.__private_update(dico)
else:
raise ArgumentTypeError(
'NLEInfos requires arguments of type dict or a list of param=value, \
see help(NLEInfos)!')
self.__private_update(kwargs)
def __private_update(self, dictionary):
"""
Doc de __private_update
"""
for key, val in dictionary.items():
if key in self.sol:
if isinstance(self.sol[key], val.__class__):
self.sol[key] = val
else:
raise ArgumentTypeError('NLEInfos ' + key + ' must be of type: ' + \
str(self.sol[key].__class__) + \
' instead of type: ' + str(val.__class__) + '!')
else:
self.sol[key] = val
def get(self, key):
"""
Doc de get
"""
try:
value = self.sol[key]
except KeyError as err:
raise OptionNameError(key + ' is not a NLEInfos, see help(NLEInfos)!\n')
return value
def __str__(self):
"""
Doc de __str__
"""
keys = sorted(self.sol.keys())
string = self.__class__.__name__ + ':\n\n'
for key in keys:
string = string + "\t" + key.ljust(20) + ' = ' + str(self.sol[key]) + '\n'
return string
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]class Options(options.Options):
"""
Options for :meth:`nutopy.nle.solve`.
Attributes
----------
Display : :class:`str`, {'on', 'off'}, default is ``on``
Display iterations and results or not
MaxFEval : :class:`int`, default is ``2000``
Maximum number of calls to the function ``f``
SolverMethod : :class:`str`, default is ``'hybrj'``
Solver method
- From nutopy package [1]_:
- `hybrj`
- From :func:`scipy.optimize.root`:
- `hybr` : this method is equivalent to `hybrj`.
- `lm`
- `broyden1`
- `broyden2`
- `anderson`
- `linearmixing`
- `diagbroyden`
- `excitingmixing`
- `krylov`
- `df-sane`
ScipyRootOptions : :class:`dict`
A dictionary of solver options for python scipy solver.
E.g. `xtol` or `maxiter`, see show_options() for details:
.. highlight:: python
>>> from scipy import optimize
>>> print(optimize.show_options(solver='root'))
TolX : :class:`float`, default is ``1e-8``
Relative tolerance between iterates.
Examples
--------
>>> from nutopy import nle
Constructor usage
>>> options = nle.Options()
>>> options = nle.Options(Display='off', TolX=1e-8)
>>> options = nle.Options({'Display' : 'off', 'TolX' : 1e-8})
Update
>>> options.update(Display='on')
Get
>>> solver = options.get('Display')
References
----------
.. [1] M. J. D. Powell, A Hybrid Method for Nonlinear Equations.
Numerical Methods for Nonlinear Algebraic Equations,
P. Rabinowitz, editor. Gordon and Breach, 1970.
"""
def __init__(self, *args, **kwargs):
"""
Doc de __init__
Parameters
----------
args : tuple
Anonymous arguments
kwargs : dictionary
Named arguments
"""
self.name = 'nle.solve'
# On donne les valeurs par defaut de options
self.options = dict(Display='on',
MaxFEval=2000,
SolverMethod='hybrj',
ScipyRootOptions=dict(),
TolX=1e-8)
self.update(*args, **kwargs)
def getAllOptions(self):
"""
Doc de getAllOptions
"""
# Cette methode permet de regrouper les options a donner en parametre aux
# fonctions fortran faisant l interface
dw = [self.options['TolX']]
iw = [self.options['MaxFEval']]
lsw = [len(self.options['Display']),
len(self.options['SolverMethod'])]
sw = self.options['Display'] + \
self.options['SolverMethod']
return [dw, iw, lsw, sw]