"""
path module description
"""
import numpy as np
from nutopy.mod_pathsolve import *
from nutopy.errors import *
from nutopy import options
from nutopy.tools import *
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
def finitediffpath(fun, x, l, *args, **kwargs):
veps = np.finfo(float).eps
n = x.size
m = l.size
j = np.zeros((n,n+m))
f = fun(x, l, *args, **kwargs)
for i in range(0,n):
t = np.sqrt(veps)*np.sqrt(np.maximum(1.0, np.abs(x[i])))
xi = x[i]
x[i] = xi + t
j[:,i] = (fun(x, l, *args, **kwargs)-f)/t
x[i] = xi
for i in range(0,m):
t = np.sqrt(veps)*np.sqrt(np.maximum(1.0, np.abs(l[i])))
li = l[i]
l[i] = li + t
j[:,n+i]= (fun(x, l, *args, **kwargs)-f)/t
l[i] = li
return j
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def solve(f, x0, pars0, parsf, *, df=None, options=None, callback=None, \
args=(), kwargs={}):
"""
Solver of one-parameter family of nonlinear equations:
.. math::
f(x, \lambda) = 0, \quad \lambda \in [\lambda_0, \lambda_f]
with :math:`f : \mathbf{R}^n \\times \mathbf{R} \\rightarrow \mathbf{R}^n`.
The solver needs an initial pair :math:`(x_0, \lambda_0)` such that
:math:`f(x_0, \lambda_0)=0`, that is the pair is a ``zero`` of :math:`f`.
The goal is to compute a path of zeros :math:`(x_n, \lambda_n)` such that the
final pair denoted :math:`(x_N, \lambda_N)` satisfies :math:`f(x_N, \lambda_N)=0`
and :math:`\lambda_N=\lambda_f`. The final point :math:`x_N` is then denoted
:math:`x_f`.
The argument :math:`\lambda` is refered as a parameter.
The solver admits a vector of parameters denoted ``pars`` instead of a scalar
parameter. In that case, the solver solves
.. math::
f(x, (1-\lambda)\, \\texttt{pars}_0 + \lambda\, \\texttt{pars}_f) = 0,
\quad \lambda \in [0, 1],
and the goal is to compute a path of zeros whose extremities satisfy
:math:`f(x_0, \\texttt{pars}_0) = 0` and :math:`f(x_f, \\texttt{pars}_f) = 0`.
Notes
-----
The differential path following algorithm of the :meth:`nutopy.path.solve` method is
partially described in Ref [1]_. The method is a Predictor-Corrector method [2]_
with high-order Runge-Kutta schemes for the prediction and a simplified
Newton for the correction. Only the prediction is descibred in [1]_.
Parameters
----------
f : callable
A vector function representing ``f(x, pars)``.
x0 : float vector or :class:`float`
Initial solution, that is ``f(x0, pars0) = 0``.
pars0 : float vector or :class:`float`
Initial parameters.
parsf : float vector or :class:`float`
Final parameters, the goal is to compute xf such that ``F(xf, parsf) = 0``.
df : callable, optional, key only
The Jacobian ``f'(x, pars)`` which returns a tuple of float matrices.
It must return a tuple ``(dfdx, dfdpars)`` where ``dfdx = df/dx(x, pars)``
and ``dfdpars = df/dpars(x,pars)``.
``dfdx`` has shape ``(n,n)`` where ``n`` is the size of ``x`` and
``dfdpars`` has shape ``(n,m)`` where ``m`` is the size of ``pars``.
options : :meth:`nutopy.path.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 ``path.solve`` method stops if ``status < 0`` at some iteration.
The argument ``infos`` is a structure containing:
- ``infos.x`` : ``x`` at current iteration;
- ``infos.pars`` : ``pars`` at current iteration;
- ``infos.f`` : ``f(x,pars)`` at current iteration.
**Advice.**
Return a status lower than -10 to avoid any confusion with status below.
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.xf`` or ``sol.get('xf')`` (float vector) : final solution ``xf``;
- ``sol.parsf`` (float vector) : final parameters ``parsf``;
- ``sol.xout`` (float matrix) : solutions along the path of zeros,
ie ``f(xout[k], parsout[k]) = 0``;
- ``sol.parsout`` (float matrix) : parameters at each computed step;
- ``sol.fnorm`` (float vector) : norm of the function ``f`` at each computed step;
- ``sol.status`` (:class:`int`) : termination status of the solver
- if ``status== 0`` then `MaxSf` is too small;
- if ``status==-1`` then input is not consistent;
- if ``status==-2`` then larger `MaxSteps` is needed;
- if ``status==-3`` then step size became too small;
- if ``status==-4`` then problem is problably stiff;
- if ``status==-5`` then ``|f|`` is greater than `NormeSfunMax`;
- if ``status==-6`` then the homotopic parameter do not
evolve relatively w.r.t. `TolHomParEvolution`;
- if ``status==-7`` then a turning point occurs.
- ``sol.success`` (:class:`bool`) : whether or not the solver exited successfully;
- ``sol.message`` (:class:`str`) : description of the cause of the termination;
- ``sol.sout`` (float vector) : arc length at each computed step;
- ``sol.vout`` (float matrix) : ``vout[:,k] = c'[:,k]`` where ``c = (x, pars)``;
- ``sol.dets`` (float vector) : ``det(f'(c(s)),c'(s))`` at each computed step;
- ``sol.ps`` (float vector) :``( c'(s_old) | c'(s) )`` at each computed step.
See Also
--------
nutopy.path.Options
nutopy.nle.solve
nutopy.errors
Notes
-----
- Use :class:`numpy.ndarray` to define float vectors.
- You can use a Newton solver (see :meth:`nutopy.nle.solve`) to get ``x0``.
- You can define a unique function ``f`` coding f and its derivative:
.. highlight:: python
>>> from nutopy.tools import *
>>> def df(x, dx, p, dp, a):
... dy = np.zeros(x.size)
... dy[0] = a * dx[0] - dp[0]
... dy[1] = np.exp(x[0] + x[1]) * ( dx[0] + dx[1] ) - dp[1]
... return dy
>>> @tensorize(df, tvars=(1, 2))
... def f(x, p, a): # p is pars
... y = np.zeros(x.size)
... y[0] = a * x[0] - p[0]
... y[1] = np.exp(x[0] + x[1]) - p[1]
... return y
Then, call ``path.solve`` with ``df=f``:
.. highlight:: python
>>> from nutopy import path
>>> x0 = np.array([0.0, 0.0])
>>> p0 = np.array([0.0, 1.0])
>>> pf = np.array([1.0, 2.0])
>>> a = 2.0
>>> sol = path.solve(f, x0, p0, pf, args=a, df=f)
>>> # or equivalently: df is not required if f is tensorized
>>> sol = path.solve(f, x0, p0, pf, args=a)
Examples
--------
>>> import numpy as np
>>> from nutopy import path
Example 1: scalar
>>> def f(x, l):
... return 2.0*x-l
>>> x0 = 0.0
>>> l0 = 0.0
>>> lf = 1.0
>>> sol = path.solve(f, x0, l0, lf)
Example 2: scalar, with numpy
>>> def f(x, l):
... y = np.zeros(x.size)
... y[0] = 2.0*x[0] - l[0]
... return y
>>> x0 = np.array([0.0])
>>> l0 = np.array([0.0])
>>> lf = np.array([1.0])
>>> sol = path.solve(f, x0, l0, lf)
Example 3: with Jacobian, options and optional argument
>>> def f(x, l, a):
... y = np.zeros(x.size)
... y[0] = a*x[0]-l[0]
... y[1] = np.exp(x[0]+x[1])-1.0
... return y
>>> def df(x, l, a):
... n = x.size
... m = l.size
... dfdx = np.zeros((n,n))
... dfdl = np.zeros((n,m))
...
... dfdx[0,0] = a
... dfdx[1,0] = np.exp(x[0]+x[1])
... dfdx[1,1] = np.exp(x[0]+x[1])
...
... dfdl[0,0] = -1.0
... return (dfdx, dfdl)
>>> opt = path.Options(ODESolver='dopri5', MaxIterCorrection=0);
>>> a = 2.0
>>> x0 = np.array([0.0, 0.0])
>>> l0 = 0.0
>>> lf = 1.0
>>> sol = path.solve(f, x0, l0, lf, options=opt, args=a, df=df)
Example 4: homotopy on two parameters with Jacobian
>>> def f(x, l, a):
... y = np.zeros(x.size)
... y[0] = a * x[0] - l[0]
... y[1] = np.exp(x[0] + x[1]) - l[1]
... return y
>>> def df(x, l, a):
... n = x.size
... m = l.size
... dfdx = np.zeros((n, n))
... dfdl = np.zeros((n, m))
...
... dfdx[0, 0] = a
... dfdx[1, 0] = np.exp(x[0] + x[1])
... dfdx[1, 1] = np.exp(x[0] + x[1])
...
... dfdl[0, 0] = -1.0
... dfdl[1, 1] = -1.0
... return (dfdx, dfdl)
>>> x0 = np.array([0.0, 0.0])
>>> l0 = np.array([0.0, 1.0])
>>> lf = np.array([1.0, 2.0])
>>> a = 2.0
>>> sol = path.solve(f, x0, l0, lf, args=a, df=df)
Example 5: use callback to stop homotopy
>>> def f(x, l):
... y = np.zeros(x.size)
... y[0] = 2.0 * x[0] - l[0]
... return y
>>> def mycallback(infos):
... if(infos.x[0]>=1.0/4.0): # we stop if x >= 1/4
... return -11
... return 0
>>> x0 = np.array([0.0])
>>> l0 = np.array([0.0])
>>> lf = np.array([1.0])
>>> sol = path.solve(f, x0, l0, lf, callback=mycallback)
>>> print(sol.status)
-11
References
----------
.. [1] J.-B. Caillau, O. Cots and J. Gergaud,
Differential continuation for regular optimal control problems,
Optimization Methods and Software, 27 (2011), no. 2, 177–196.
.. [2] E. Allgower & K. Georg, Introduction to numerical continuation methods,
vol. 45 of Classics in Applied Mathematics, Soc. for Industrial and Applied Math.,
Philadelphia, PA, USA, (2003), xxvi+388.
Raises
------
ArgumentTypeError
Variable x0 must be a one dimensional numpy array or a float.
ArgumentTypeError
Variable pars0 must be a one dimensional numpy array or a float.
ArgumentTypeError
Variable parsf must be a one dimensional numpy array or a float.
ArgumentTypeError
Variable options must be a Options object.
ArgumentDimensionError
Variable x0 must be a one dimensional numpy array or a float.
ArgumentDimensionError
Variable pars0 must be a one dimensional numpy array or a float.
ArgumentDimensionError
Variable parsf must be a one dimensional numpy array or a float.
ArgumentDimensionError
Variable parsf must have the same size than pars0.
"""
# -----------------------------------------------------------------------------------
#
# Cast of input variables
#
# -----------------------------------------------------------------------------------
if(not isinstance(x0, float)):
x0 = np.asarray(x0)
if(not isinstance(pars0, float)):
pars0 = np.asarray(pars0)
if(not isinstance(parsf, float)):
parsf = np.asarray(parsf)
# -----------------------------------------------------------------------------------
#
# Tests on input variables
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
if options is None:
options = Options()
# -----------------------------------------------------------------------------------
#Variable x0 must be a vector
# -----------------------------------------------------------------------------------
# Variable x0 must be a vector
msgErreur = '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(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0,float)):
n = 1
else :
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
msgErreur = 'Variable pars0 must be a one dimensional numpy array or a float.';
if(isinstance(pars0,np.ndarray)):
#Variable pars0 must be a vector
if(pars0.ndim==1):
np0 = pars0.size
elif(pars0.ndim==2):
(lig,col) = pars0.shape
if(lig==1 or col==1):
np0 = pars0.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(pars0,float)):
np0 = 1
else :
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
msgErreur = 'Variable parsf must be a one dimensional numpy array or a float.';
if(isinstance(parsf,np.ndarray)):
#Variable parsf must be a vector
if(parsf.ndim==1):
npf = parsf.size
elif(parsf.ndim==2):
(lig,col) = parsf.shape
if(lig==1 or col==1):
npf = parsf.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(parsf,float)):
npf = 1
else :
raise ArgumentTypeError(msgErreur)
#Test if np0 = npf
if(np0!=npf):
raise ArgumentDimensionError('Variable parsf must have the same size than pars0.')
# -----------------------------------------------------------------------------------
#Variable options must be a Options object
msgErreur = 'Variable options must be a Options object!';
if(not isinstance(options,Options)):
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#
# create callback function passed to fortran routine
#
display = options.get('Display')=='on'
disp_iter = options.get('DispIter')
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, 2), x0, pars0, *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: finite diff
myjac = lambda x, pars: finitediffpath(f, x, pars, *args, **kwargs)
elif jac_case == 2: # f is tensorized
def myjac(x, pars):
n = x.size
m = pars.size
dfdx = np.zeros((n,n))
dfdp = np.zeros((n,m))
Idx = np.eye(n)
Idp = np.eye(m)
zx = np.zeros(n)
zp = np.zeros(m)
for i in range(0, n): _, dfdx[:,i] = df((x, Idx[:, i]), (pars, zp), *args, **kwargs)
for i in range(0, m): _, dfdp[:,i] = df((x, zx), (pars, Idp[:, i]), *args, **kwargs)
return dfdx, dfdp
else: # jac_case = 1: the Jacobian is given
myjac = lambda x, pars: df(x, pars, *args, **kwargs)
myfun = lambda x, pars: f(x, pars, *args, **kwargs)
if callback is None:
def callbackPATH(x, pars, fx, hompar, arclength, dets, dots, Niter):
if (display) and (disp_iter>=0) and ((Niter%disp_iter==0) or (Niter==1)):
callbackPrint(x, fx, hompar, arclength, dets, dots, Niter)
return 0
else:
def callbackPATH(x, pars, fx, hompar, arclength, dets, dots, Niter):
if (display) and (disp_iter>=0) and ((Niter%disp_iter==0) or (Niter==1)):
callbackPrint(x, fx, hompar, arclength, dets, dots, Niter)
infos_dict = dict(x=x, pars=pars, f=fx)
infos = PATHInfos(infos_dict)
if callback.__code__.co_argcount > 1:
return callback(infos, *args, **kwargs)
else:
return callback(infos)
# -----------------------------------------------------------------------------------
# on passe du tuple (df/dx, df/dpars) a un array
def myjacAux(x, pars):
n = x.size
m = pars.size
j = np.zeros((n,n+m))
(dfdx, dfdpars) = myjac(x, pars)
j[:,0:n] = dfdx
j[:,n:] = dfdpars
return j
# -----------------------------------------------------------------------------------
if df is None:
jac_final = myjac
else:
jac_final = myjacAux
# -----------------------------------------------------------------------------------
#
if(not isinstance(x0,float)):
x0.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
#
# prepare call to Fortran solver
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[dw, iw, lsw, sw] = options.getAllOptions()
ndw = len(dw)
niw = len(iw)
nsw = len(lsw)
nswint = len(sw)
swint = [ord(c) for c in sw] #Code ascii de la chaine de caractere
# parspan
nparspan = np0
mparspan = 2
parspan = np.zeros((nparspan,mparspan),dtype='d',order='F')
parspan[:,0] = pars0
parspan[:,1] = parsf
# -----------------------------------------------------------------------------------
[dimepath,flag,messageFortran] = mod_pathsolve.pathsolvepy(myfun,jac_final,x0, \
parspan,callbackPATH,dw,iw,swint,lsw,n,nparspan,mparspan,ndw,niw,nsw,nswint)
# -----------------------------------------------------------------------------------
#
# post computations, prepare outputs
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
parsout_for = mod_pathsolve.parout
xout_for = mod_pathsolve.xout
sout = mod_pathsolve.sout
viout = mod_pathsolve.viout # viout.shape = (n+1, len(sout))
dets = mod_pathsolve.dets
fnorm = mod_pathsolve.normf
ps = mod_pathsolve.ps
# -----------------------------------------------------------------------------------
# on tensorise : parsout, xout, viout
if(not isinstance(x0,float)):
x0shape = x0.shape
xout = np.reshape(xout_for[:, 0], x0shape)[None]
else:
xout = xout_for[0, 0][None]
for i in range(dimepath-1):
if(not isinstance(x0,float)):
x = np.reshape(xout_for[:, i+1], x0shape)
else:
x = xout_for[0, i+1]
xout = np.concatenate((xout, x[None]), axis=0)
if(not isinstance(pars0,float)):
pars0shape = pars0.shape
parsout = np.reshape(parsout_for[:, 0], pars0shape)[None]
else:
parsout = parsout_for[0, 0][None]
for i in range(dimepath-1):
if(not isinstance(pars0,float)):
p = np.reshape(parsout_for[:, i+1], pars0shape)
else:
p = parsout_for[0, i+1]
parsout = np.concatenate((parsout, p[None]), axis=0)
# -----------------------------------------------------------------------------------
parsfin = parsout[-1]
xf = xout[-1]
# -----------------------------------------------------------------------------------
message = messageFortran.strip().decode()
status = flag
if(status==1):
success=True
else:
success=False
# -----------------------------------------------------------------------------------
sol = dict(xf=xf, parsf=parsfin, parsout=parsout, xout=xout, sout=sout, vout=viout, \
dets=dets, fnorm=fnorm, ps=ps, status=status, success=success, message=message)
pathsol = PATHSol(sol)
# -----------------------------------------------------------------------------------
# display outputs
if display:
print('\n Results of the path solver method:\n')
print(' xf = ', pathsol.xf);
print(' parsf = ', pathsol.parsf);
print(' |f(xf,parsf)| = ', pathsol.fnorm[-1]);
print(' steps = ', pathsol.sout.size);
print(' status = ', pathsol.status);
print(' success = ', pathsol.success, '\n');
print(' ' + pathsol.message + '\n')
# -----------------------------------------------------------------------------------
# reset counter
callbackPrint.counter=0
return pathsol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
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, fx, hompar, arclength, dets, dots, Niter):
if(callbackPrint.counter==0):
# do not modify the following line
print('\n Calls |f(x,pars)| |x| Homotopic param \
Arclength s det A(s) Dot product \n ')
print('{0:10}'.format(Niter) + '{0:16.8e}'.format(np.linalg.norm(fx)) + \
'{0:19.11e}'.format(np.linalg.norm(x)) + '{0:19.11e}'.format(hompar) \
+ '{0:16.8e}'.format(arclength) + '{0:16.8e}'.format(dets) + \
'{0:16.8e}'.format(dots))
callbackPrint.counter += 1
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class PATHSol:
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
self.sol = dict()
self.update(*args, **kwargs)
self.xf = self.sol['xf']
self.parsf = self.sol['parsf']
self.parsout = self.sol['parsout']
self.xout = self.sol['xout']
self.sout = self.sol['sout']
self.vout = self.sol['vout']
self.dets = self.sol['dets']
self.fnorm = self.sol['fnorm']
self.ps = self.sol['ps']
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('PATHSol requires arguments of type dict or \
a list of param=value, see help(PATHSol)!')
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('PATHSol ' + 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 PATHSol, see help(PATHSol)!\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 PATHInfos:
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
self.sol = dict()
self.update(*args, **kwargs)
self.x = self.sol['x']
self.pars = self.sol['pars']
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('PATHInfos requires arguments of type dict or \
a list of param=value, see help(PATHInfos)!')
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('PATHInfos ' + 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 PATHInfos, see help(PATHInfos)!\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.path.solve`.
Attributes
----------
DispIter : :class:`int`, default is ``1``
Display at iteration k such that mod(k,DispIter)==0
Display : :class:`str`, {'on', 'off'}, default is ``on``
Display iterations and results or not
DoSavePath : :class:`str`, {'off', 'on'}, default is ``off``
Save path during the homotopy in a file "SavedPath(n).py"
MaxIterCorrection : :class:`int`, default is ``7``
Maximum number of iterations during correction
MaxArcLength : :class:`float`, default is ``1e5``
Maximum arc length
MaxFunNorm : :class:`float`, default is ``1e-1``
Maximal norm of the function F during homotopy
MaxSteps : :class:`int`, default is ``10000``
Maximum number of homotopy steps
MaxStepSize : :class:`float`, default is ``0``
Maximum step size during homotopy.
If it is equal to 0 then the integrator has no constraint on the maximal step size.
MaxStepSizeHomPar : :class:`float`, default is ``0``
Maximum step size for the homotopic parameter.
If 0 then no constraint on the maximal step size.
ODESolver : :class:`str`, default is ``dopri5``
Integrator name for homotopy. See Refs. [1]_ and [2]_.
The integrators are interface of Fortran softwares available
`here <http://www.unige.ch/~hairer/software.html>`_.
- Explicit:
dopri5, dop853
- Implicit:
radau5, radau9, radau13, radau (adaptative order 5, 9 and 13)
StopAtTurningPoint : :class:`int`, {0, 1}, default is ``0``
Stop or not after a turning point.
TolOdeAbs : :class:`float`, default is ``1e-10``
Absolute error tolerance.
TolOdeRel : :class:`float`, default is ``1e-8``
Relative error tolerance.
TolHomparfinalStep : :class:`float`, default is ``1e-8``
Absolute Dense output tolerance.
Absolute tolerance to detect if the final homotopic parameter has been reached.
TolHomParEvolution : :class:`float`, default is ``1e-8``
Relative Dense output tolerance.
The homotopy stops when the homotopic parameter do not
evolve relatively, iteration per iteration, during ten following steps.
TolXCorrectionStep : :class:`float`, default is ``1e-8``
Relative tolerance during correction.
Examples
--------
>>> from nutopy import path
Constructor usage
>>> options = path.Options()
>>> options = path.Options(Display='off', MaxSf=100)
>>> options = path.Options({'Display' : 'off', 'MaxSf' : 100})
Update
>>> options.update(Display='on')
Get
>>> solver = options.get('Display')
References
----------
.. [1] E. Hairer, S. P. Norsett & G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems,
vol 8 of Springer Serie in Computational Mathematics, Springer-Verlag, second edn (1993).
.. [2] E. Hairer & G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems,
vol 14 of Springer Serie in Computational Mathematics, Springer-Verlag, second edn (1996).
"""
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
self.name = 'path.solve'
#On donne les valeurs par defaut de options
self.options = dict(DispIter = 1,
Display = 'on',
DoSavePath = 'off',
MaxArcLength = 1e5,
MaxFunNorm = 1e-1,
MaxSteps = 10000,
MaxStepSize = 0.0,
ODESolver = 'dopri5',
StopAtTurningPoint = 0,
TolOdeAbs = 1e-10,
TolOdeRel = 1e-8,
TolHomparfinalStep = 1e-8,
TolHomParEvolution = 1e-8,
TolXCorrectionStep = 1e-8,
MaxIterCorrection = 7,
MaxStepSizeHomPar = 0.0)
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
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['MaxArcLength'] ,
self.options['MaxFunNorm'] ,
self.options['MaxStepSize'] ,
self.options['MaxStepSizeHomPar'] ,
self.options['TolOdeAbs'] ,
self.options['TolOdeRel'] ,
self.options['TolHomparfinalStep'] ,
self.options['TolHomParEvolution'] ,
self.options['TolXCorrectionStep'] ]
iw = [ self.options['DispIter'] ,
self.options['MaxIterCorrection'] ,
self.options['MaxSteps'] ,
self.options['StopAtTurningPoint'] ]
lsw = [ len(self.options['Display']) ,
len(self.options['DoSavePath']) ,
len(self.options['ODESolver']) ]
sw = self.options['Display'] + \
self.options['DoSavePath'] + \
self.options['ODESolver']
return [dw, iw, lsw, sw]