from __future__ import absolute_import
import numpy as nm
import numpy.linalg as nla
from sfepy.base.base import output, get_default, pause, Struct
from sfepy.base.log import Log, get_logging_conf
from sfepy.base.timing import Timer
from sfepy.solvers.solvers import OptimizationSolver
import scipy.optimize as sopt
import scipy.optimize.linesearch as linesearch
import six
from six.moves import range
[docs]
def conv_test(conf, it, of, of0, ofg_norm=None):
"""
Returns
-------
flag : int
* -1 ... continue
* 0 ... small OF -> stop
* 1 ... i_max reached -> stop
* 2 ... small OFG -> stop
* 3 ... small relative decrase of OF
"""
status = -1
output('opt: iter: %d, of: %e (||ofg||: %e)' % (it, of, ofg_norm))
if (abs(of) < conf.eps_of):
status = 0
elif ofg_norm and (ofg_norm < conf.eps_ofg):
status = 2
elif (it > 0) and (abs(of0 - of) < (conf.eps_rd * abs(of0))):
status = 3
if (status == -1) and (it >= conf.i_max):
status = 1
return status
[docs]
def wrap_function(function, args):
ncalls = [0]
times = []
timer = Timer()
def function_wrapper(x):
ncalls[0] += 1
timer.start()
out = function(x, *args)
times.append(timer.stop())
return out
return ncalls, times, function_wrapper
[docs]
def check_gradient(xit, aofg, fn_of, delta, check):
dofg = nm.zeros_like(aofg)
xd = xit.copy()
for ii in range(xit.shape[0]):
xd[ii] = xit[ii] + delta
ofp = fn_of(xd)
xd[ii] = xit[ii] - delta
ofm = fn_of(xd)
xd[ii] = xit[ii]
dofg[ii] = 0.5 * (ofp - ofm) / delta
output('**********', ii, aofg[ii], dofg[ii])
diff = abs(aofg - dofg)
aux = nm.concatenate((aofg[:,nm.newaxis], dofg[:,nm.newaxis],
diff[:,nm.newaxis]), 1)
output(aux)
output(nla.norm(diff, nm.Inf))
aofg.tofile('aofg.txt', ' ')
dofg.tofile('dofg.txt', ' ')
diff.tofile('diff.txt', ' ')
if check == 2:
import pylab
pylab.plot(aofg)
pylab.plot(dofg)
pylab.legend(('analytical', 'finite difference'))
pylab.show()
pause('gradient checking done')
[docs]
class FMinSteepestDescent(OptimizationSolver):
"""
Steepest descent optimization solver.
"""
name = 'opt.fmin_sd'
_parameters = [
('i_max', 'int', 10, False,
'The maximum number of iterations.'),
('eps_rd', 'float', 1e-5, False,
'The relative delta of the objective function.'),
('eps_of', 'float', 1e-4, False,
'The tolerance for the objective function.'),
('eps_ofg', 'float', 1e-8, False,
'The tolerance for the objective function gradient.'),
('norm', 'numpy norm', nm.Inf, False,
'The norm to be used.'),
('ls', 'bool', True, False,
'If True, use a line-search.'),
('ls_method', "{'backtracking', 'full'}", 'backtracking', False,
'The line-search method.'),
('ls_on', 'float', 0.99999, False,
"""Start the backtracking line-search by reducing the step, if
:math:`||f(x^i)|| / ||f(x^{i-1})||` is larger than `ls_on`."""),
('ls0', '0.0 < float < 1.0', 1.0, False,
'The initial step.'),
('ls_red', '0.0 < float < 1.0', 0.5, False,
'The step reduction factor in case of correct residual assembling.'),
('ls_red_warp', '0.0 < float < 1.0', 0.1, False,
"""The step reduction factor in case of failed residual assembling
(e.g. the "warp violation" error caused by a negative volume
element resulting from too large deformations)."""),
('ls_min', '0.0 < float < 1.0', 1e-5, False,
'The minimum step reduction factor.'),
('check', '0, 1 or 2', 0, False,
"""If >= 1, check the tangent matrix using finite differences. If 2,
plot the resulting sparsity patterns."""),
('delta', 'float', 1e-6, False,
r"""If `check >= 1`, the finite difference matrix is taken as
:math:`A_{ij} = \frac{f_i(x_j + \delta) - f_i(x_j - \delta)}{2
\delta}`."""),
('output', 'function', None, False,
"""If given, use it instead of :func:`output()
<sfepy.base.base.output()>` function."""),
('yscales', 'list of str', ['linear', 'log', 'log', 'linear'], False,
'The list of four convergence log subplot scales.'),
('log', 'dict or None', None, False,
"""If not None, log the convergence according to the configuration in
the following form: ``{'text' : 'log.txt', 'plot' : 'log.pdf'}``.
Each of the dict items can be None."""),
]
def __init__(self, conf, **kwargs):
OptimizationSolver.__init__(self, conf, **kwargs)
conf = self.conf
log = get_logging_conf(conf)
conf.log = log = Struct(name='log_conf', **log)
conf.is_any_log = (log.text is not None) or (log.plot is not None)
if conf.is_any_log:
self.log = Log([[r'$||\Psi||$'], [r'$||\nabla \Psi||$'],
[r'$\alpha$'], ['iteration']],
xlabels=['', '', 'all iterations', 'all iterations'],
yscales=conf.yscales,
is_plot=conf.log.plot is not None,
log_filename=conf.log.text,
formats=[['%.8e'], ['%.3e'], ['%.3e'], ['%d']])
else:
self.log = None
def __call__(self, x0, conf=None, obj_fun=None, obj_fun_grad=None,
status=None, obj_args=None):
conf = get_default(conf, self.conf)
obj_fun = get_default(obj_fun, self.obj_fun)
obj_fun_grad = get_default(obj_fun_grad, self.obj_fun_grad)
status = get_default(status, self.status)
obj_args = get_default(obj_args, self.obj_args)
if conf.output:
globals()['output'] = conf.output
output('entering optimization loop...')
nc_of, tt_of, fn_of = wrap_function(obj_fun, obj_args)
nc_ofg, tt_ofg, fn_ofg = wrap_function(obj_fun_grad, obj_args)
timer = Timer()
time_stats = {'of' : tt_of, 'ofg': tt_ofg, 'check' : []}
ofg = None
it = 0
xit = x0.copy()
while 1:
of = fn_of(xit)
if it == 0:
of0 = ofit0 = of_prev = of
of_prev_prev = of + 5000.0
if ofg is None:
ofg = fn_ofg(xit)
if conf.check:
timer.start()
check_gradient(xit, ofg, fn_of, conf.delta, conf.check)
time_stats['check'].append(timer.stop())
ofg_norm = nla.norm(ofg, conf.norm)
ret = conv_test(conf, it, of, ofit0, ofg_norm)
if ret >= 0:
break
ofit0 = of
##
# Backtrack (on errors).
alpha = conf.ls0
can_ls = True
while 1:
xit2 = xit - alpha * ofg
aux = fn_of(xit2)
if self.log is not None:
self.log(of, ofg_norm, alpha, it)
if aux is None:
alpha *= conf.ls_red_warp
can_ls = False
output('warp: reducing step (%f)' % alpha)
elif conf.ls and conf.ls_method == 'backtracking':
if aux < of * conf.ls_on: break
alpha *= conf.ls_red
output('backtracking: reducing step (%f)' % alpha)
else:
of_prev_prev = of_prev
of_prev = aux
break
if alpha < conf.ls_min:
if aux is None:
raise RuntimeError('giving up...')
output('linesearch failed, continuing anyway')
break
# These values are modified by the line search, even if it fails
of_prev_bak = of_prev
of_prev_prev_bak = of_prev_prev
if conf.ls and can_ls and conf.ls_method == 'full':
output('full linesearch...')
alpha, fc, gc, of_prev, of_prev_prev, ofg1 = \
linesearch.line_search(fn_of,fn_ofg,xit,
-ofg,ofg,of_prev,of_prev_prev,
c2=0.4)
if alpha is None: # line search failed -- use different one.
alpha, fc, gc, of_prev, of_prev_prev, ofg1 = \
sopt.line_search(fn_of,fn_ofg,xit,
-ofg,ofg,of_prev_bak,
of_prev_prev_bak)
if alpha is None or alpha == 0:
# This line search also failed to find a better
# solution.
ret = 3
break
output(' -> alpha: %.8e' % alpha)
else:
if conf.ls_method == 'full':
output('full linesearch off (%s and %s)'
% (conf.ls, can_ls))
ofg1 = None
if self.log is not None:
self.log.plot_vlines(color='g', linewidth=0.5)
xit = xit - alpha * ofg
if ofg1 is None:
ofg = None
else:
ofg = ofg1.copy()
for key, val in six.iteritems(time_stats):
if len(val):
output('%10s: %7.2f [s]' % (key, val[-1]))
it = it + 1
output('status: %d' % ret)
output('initial value: %.8e' % of0)
output('current value: %.8e' % of)
output('iterations: %d' % it)
output('function evaluations: %d in %.2f [s]'
% (nc_of[0], nm.sum(time_stats['of'])))
output('gradient evaluations: %d in %.2f [s]'
% (nc_ofg[0], nm.sum(time_stats['ofg'])))
if self.log is not None:
self.log(of, ofg_norm, alpha, it)
if conf.log.plot is not None:
self.log(save_figure=conf.log.plot, finished=True)
else:
self.log(finished=True)
if status is not None:
status['log'] = self.log
status['status'] = status
status['of0'] = of0
status['of'] = of
status['it'] = it
status['nc_of'] = nc_of[0]
status['nc_ofg'] = nc_ofg[0]
status['time_stats'] = time_stats
return xit
[docs]
class ScipyFMinSolver(OptimizationSolver):
"""
Interface to SciPy optimization solvers scipy.optimize.fmin_*.
"""
name = 'nls.scipy_fmin_like'
_i_max_name = {
'fmin' : 'maxiter',
'fmin_bfgs' : 'maxiter',
'fmin_cg' : 'maxiter',
'fmin_cobyla' : 'maxfun',
'fmin_l_bfgs_b' : 'maxfun',
'fmin_ncg' : 'maxiter',
'fmin_powell' : 'maxiter',
'fmin_slsqp' : 'iter',
'fmin_tnc' : 'maxfun',
}
_has_grad = ('fmin_bfgs', 'fmin_cg', 'fmin_l_bfgs_b', 'fmin_ncg',
'fmin_slsqp', 'fmin_tnc')
_parameters = [
('method',
'{%s}' % ', '.join(sorted(repr(ii) for ii in _i_max_name.keys())),
'fmin', False,
'The actual optimization method to use.'),
('i_max', 'int', 10, False,
'The maximum number of iterations.'),
('*', '*', None, False,
'Additional parameters supported by the method.'),
]
def __init__(self, conf, **kwargs):
OptimizationSolver.__init__(self, conf, **kwargs)
self.set_method(self.conf)
[docs]
def set_method(self, conf):
import scipy.optimize as so
try:
solver = getattr(so, conf.method)
except AttributeError:
raise ValueError('scipy solver %s does not exist!' % conf.method)
self.solver = solver
def __call__(self, x0, conf=None, obj_fun=None, obj_fun_grad=None,
status=None, obj_args=None):
import inspect
if conf is not None:
self.set_method(conf)
else:
conf = self.conf
obj_fun = get_default(obj_fun, self.obj_fun)
obj_fun_grad = get_default(obj_fun_grad, self.obj_fun_grad)
status = get_default(status, self.status)
obj_args = get_default(obj_args, self.obj_args)
timer = Timer(start=True)
kwargs = {self._i_max_name[conf.method] : conf.i_max,
'args' : obj_args}
if conf.method in self._has_grad:
kwargs['fprime'] = obj_fun_grad
if 'disp' in inspect.getargspec(self.solver)[0]:
kwargs['disp'] = conf.verbose
kwargs.update(self.build_solver_kwargs(conf))
out = self.solver(obj_fun, x0, **kwargs)
if status is not None:
status['time_stats'] = timer.stop()
return out