src/daComposant/daAlgorithms/Atoms/lbfgsb19hlt.py

   1 # Modification de la version 1.9.1
   2 """
   3 Functions
   4 ---------
   5 .. autosummary::
   6    :toctree: generated/
   7
   8     fmin_l_bfgs_b
   9
  10 """
  11
  12 ## License for the Python wrapper
  13 ## ==============================
  14
  15 ## Copyright (c) 2004 David M. Cooke <cookedm@physics.mcmaster.ca>
  16
  17 ## Permission is hereby granted, free of charge, to any person obtaining a
  18 ## copy of this software and associated documentation files (the "Software"),
  19 ## to deal in the Software without restriction, including without limitation
  20 ## the rights to use, copy, modify, merge, publish, distribute, sublicense,
  21 ## and/or sell copies of the Software, and to permit persons to whom the
  22 ## Software is furnished to do so, subject to the following conditions:
  23
  24 ## The above copyright notice and this permission notice shall be included in
  25 ## all copies or substantial portions of the Software.
  26
  27 ## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  28 ## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  29 ## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  30 ## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  31 ## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  32 ## FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  33 ## DEALINGS IN THE SOFTWARE.
  34
  35 ## Modifications by Travis Oliphant and Enthought, Inc. for inclusion in SciPy
  36
  37 import numpy as np
  38 from numpy import array, asarray, float64, zeros
  39 from scipy.optimize import _lbfgsb
  40 from scipy.optimize._optimize import (MemoizeJac, OptimizeResult,
  41                        _check_unknown_options, _prepare_scalar_function)
  42 from scipy.optimize._constraints import old_bound_to_new
  43
  44 from scipy.sparse.linalg import LinearOperator
  45
  46 __all__ = ['fmin_l_bfgs_b', 'LbfgsInvHessProduct']
  47
  48
  49 def fmin_l_bfgs_b(func, x0, fprime=None, args=(),
  50                   approx_grad=0,
  51                   bounds=None, m=10, factr=1e7, pgtol=1e-5,
  52                   epsilon=1e-8,
  53                   iprint=-1, maxfun=15000, maxiter=15000, disp=None,
  54                   callback=None, maxls=20):
  55     """
  56     Minimize a function func using the L-BFGS-B algorithm.
  57
  58     Parameters
  59     ----------
  60     func : callable f(x,*args)
  61         Function to minimize.
  62     x0 : ndarray
  63         Initial guess.
  64     fprime : callable fprime(x,*args), optional
  65         The gradient of `func`. If None, then `func` returns the function
  66         value and the gradient (``f, g = func(x, *args)``), unless
  67         `approx_grad` is True in which case `func` returns only ``f``.
  68     args : sequence, optional
  69         Arguments to pass to `func` and `fprime`.
  70     approx_grad : bool, optional
  71         Whether to approximate the gradient numerically (in which case
  72         `func` returns only the function value).
  73     bounds : list, optional
  74         ``(min, max)`` pairs for each element in ``x``, defining
  75         the bounds on that parameter. Use None or +-inf for one of ``min`` or
  76         ``max`` when there is no bound in that direction.
  77     m : int, optional
  78         The maximum number of variable metric corrections
  79         used to define the limited memory matrix. (The limited memory BFGS
  80         method does not store the full hessian but uses this many terms in an
  81         approximation to it.)
  82     factr : float, optional
  83         The iteration stops when
  84         ``(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps``,
  85         where ``eps`` is the machine precision, which is automatically
  86         generated by the code. Typical values for `factr` are: 1e12 for
  87         low accuracy; 1e7 for moderate accuracy; 10.0 for extremely
  88         high accuracy. See Notes for relationship to `ftol`, which is exposed
  89         (instead of `factr`) by the `scipy.optimize.minimize` interface to
  90         L-BFGS-B.
  91     pgtol : float, optional
  92         The iteration will stop when
  93         ``max{|proj g_i | i = 1, ..., n} <= pgtol``
  94         where ``pg_i`` is the i-th component of the projected gradient.
  95     epsilon : float, optional
  96         Step size used when `approx_grad` is True, for numerically
  97         calculating the gradient
  98     iprint : int, optional
  99         Controls the frequency of output. ``iprint < 0`` means no output;
 100         ``iprint = 0``    print only one line at the last iteration;
 101         ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations;
 102         ``iprint = 99``   print details of every iteration except n-vectors;
 103         ``iprint = 100``  print also the changes of active set and final x;
 104         ``iprint > 100``  print details of every iteration including x and g.
 105     disp : int, optional
 106         If zero, then no output. If a positive number, then this over-rides
 107         `iprint` (i.e., `iprint` gets the value of `disp`).
 108     maxfun : int, optional
 109         Maximum number of function evaluations. Note that this function
 110         may violate the limit because of evaluating gradients by numerical
 111         differentiation.
 112     maxiter : int, optional
 113         Maximum number of iterations.
 114     callback : callable, optional
 115         Called after each iteration, as ``callback(xk)``, where ``xk`` is the
 116         current parameter vector.
 117     maxls : int, optional
 118         Maximum number of line search steps (per iteration). Default is 20.
 119
 120     Returns
 121     -------
 122     x : array_like
 123         Estimated position of the minimum.
 124     f : float
 125         Value of `func` at the minimum.
 126     d : dict
 127         Information dictionary.
 128
 129         * d['warnflag'] is
 130
 131           - 0 if converged,
 132           - 1 if too many function evaluations or too many iterations,
 133           - 2 if stopped for another reason, given in d['task']
 134
 135         * d['grad'] is the gradient at the minimum (should be 0 ish)
 136         * d['funcalls'] is the number of function calls made.
 137         * d['nit'] is the number of iterations.
 138
 139     See also
 140     --------
 141     minimize: Interface to minimization algorithms for multivariate
 142         functions. See the 'L-BFGS-B' `method` in particular. Note that the
 143         `ftol` option is made available via that interface, while `factr` is
 144         provided via this interface, where `factr` is the factor multiplying
 145         the default machine floating-point precision to arrive at `ftol`:
 146         ``ftol = factr * numpy.finfo(float).eps``.
 147
 148     Notes
 149     -----
 150     License of L-BFGS-B (FORTRAN code):
 151
 152     The version included here (in fortran code) is 3.0
 153     (released April 25, 2011). It was written by Ciyou Zhu, Richard Byrd,
 154     and Jorge Nocedal <nocedal@ece.nwu.edu>. It carries the following
 155     condition for use:
 156
 157     This software is freely available, but we expect that all publications
 158     describing work using this software, or all commercial products using it,
 159     quote at least one of the references given below. This software is released
 160     under the BSD License.
 161
 162     References
 163     ----------
 164     * R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound
 165       Constrained Optimization, (1995), SIAM Journal on Scientific and
 166       Statistical Computing, 16, 5, pp. 1190-1208.
 167     * C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B,
 168       FORTRAN routines for large scale bound constrained optimization (1997),
 169       ACM Transactions on Mathematical Software, 23, 4, pp. 550 - 560.
 170     * J.L. Morales and J. Nocedal. L-BFGS-B: Remark on Algorithm 778: L-BFGS-B,
 171       FORTRAN routines for large scale bound constrained optimization (2011),
 172       ACM Transactions on Mathematical Software, 38, 1.
 173
 174     """
 175     # handle fprime/approx_grad
 176     if approx_grad:
 177         fun = func
 178         jac = None
 179     elif fprime is None:
 180         fun = MemoizeJac(func)
 181         jac = fun.derivative
 182     else:
 183         fun = func
 184         jac = fprime
 185
 186     # build options
 187     if disp is None:
 188         disp = iprint
 189     opts = {'disp': disp,
 190             'iprint': iprint,
 191             'maxcor': m,
 192             'ftol': factr * np.finfo(float).eps,
 193             'gtol': pgtol,
 194             'eps': epsilon,
 195             'maxfun': maxfun,
 196             'maxiter': maxiter,
 197             'callback': callback,
 198             'maxls': maxls}
 199
 200     res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
 201                            **opts)
 202     d = {'grad': res['jac'],
 203          'task': res['message'],
 204          'funcalls': res['nfev'],
 205          'nit': res['nit'],
 206          'warnflag': res['status']}
 207     f = res['fun']
 208     x = res['x']
 209
 210     return x, f, d
 211
 212
 213 def _minimize_lbfgsb(fun, x0, args=(), jac=None, bounds=None,
 214                      disp=None, maxcor=10, ftol=2.2204460492503131e-09,
 215                      gtol=1e-5, eps=1e-8, maxfun=15000, maxiter=15000,
 216                      iprint=-1, callback=None, maxls=20,
 217                      finite_diff_rel_step=None, **unknown_options):
 218     """
 219     Minimize a scalar function of one or more variables using the L-BFGS-B
 220     algorithm.
 221
 222     Options
 223     -------
 224     disp : None or int
 225         If `disp is None` (the default), then the supplied version of `iprint`
 226         is used. If `disp is not None`, then it overrides the supplied version
 227         of `iprint` with the behaviour you outlined.
 228     maxcor : int
 229         The maximum number of variable metric corrections used to
 230         define the limited memory matrix. (The limited memory BFGS
 231         method does not store the full hessian but uses this many terms
 232         in an approximation to it.)
 233     ftol : float
 234         The iteration stops when ``(f^k -
 235         f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol``.
 236     gtol : float
 237         The iteration will stop when ``max{|proj g_i | i = 1, ..., n}
 238         <= gtol`` where ``pg_i`` is the i-th component of the
 239         projected gradient.
 240     eps : float or ndarray
 241         If `jac is None` the absolute step size used for numerical
 242         approximation of the jacobian via forward differences.
 243     maxfun : int
 244         Maximum number of function evaluations.
 245     maxiter : int
 246         Maximum number of iterations.
 247     iprint : int, optional
 248         Controls the frequency of output. ``iprint < 0`` means no output;
 249         ``iprint = 0``    print only one line at the last iteration;
 250         ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations;
 251         ``iprint = 99``   print details of every iteration except n-vectors;
 252         ``iprint = 100``  print also the changes of active set and final x;
 253         ``iprint > 100``  print details of every iteration including x and g.
 254     callback : callable, optional
 255         Called after each iteration, as ``callback(xk)``, where ``xk`` is the
 256         current parameter vector.
 257     maxls : int, optional
 258         Maximum number of line search steps (per iteration). Default is 20.
 259     finite_diff_rel_step : None or array_like, optional
 260         If `jac in ['2-point', '3-point', 'cs']` the relative step size to
 261         use for numerical approximation of the jacobian. The absolute step
 262         size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
 263         possibly adjusted to fit into the bounds. For ``method='3-point'``
 264         the sign of `h` is ignored. If None (default) then step is selected
 265         automatically.
 266
 267     Notes
 268     -----
 269     The option `ftol` is exposed via the `scipy.optimize.minimize` interface,
 270     but calling `scipy.optimize.fmin_l_bfgs_b` directly exposes `factr`. The
 271     relationship between the two is ``ftol = factr * numpy.finfo(float).eps``.
 272     I.e., `factr` multiplies the default machine floating-point precision to
 273     arrive at `ftol`.
 274
 275     """
 276     _check_unknown_options(unknown_options)
 277     m = maxcor
 278     pgtol = gtol
 279     factr = ftol / np.finfo(float).eps
 280
 281     x0 = asarray(x0).ravel()
 282     n, = x0.shape
 283
 284     if bounds is None:
 285         bounds = [(None, None)] * n
 286     if len(bounds) != n:
 287         raise ValueError('length of x0 != length of bounds')
 288
 289     # unbounded variables must use None, not +-inf, for optimizer to work properly
 290     bounds = [(None if l == -np.inf else l, None if u == np.inf else u) for l, u in bounds]
 291     # LBFGSB is sent 'old-style' bounds, 'new-style' bounds are required by
 292     # approx_derivative and ScalarFunction
 293     new_bounds = old_bound_to_new(bounds)
 294
 295     # check bounds
 296     if (new_bounds[0] > new_bounds[1]).any():
 297         raise ValueError("LBFGSB - one of the lower bounds is greater than an upper bound.")
 298
 299     # initial vector must lie within the bounds. Otherwise ScalarFunction and
 300     # approx_derivative will cause problems
 301     x0 = np.clip(x0, new_bounds[0], new_bounds[1])
 302
 303     if disp is not None:
 304         if disp == 0:
 305             iprint = -1
 306         else:
 307             iprint = disp
 308
 309     sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
 310                                   bounds=new_bounds,
 311                                   finite_diff_rel_step=finite_diff_rel_step)
 312
 313     func_and_grad = sf.fun_and_grad
 314
 315     fortran_int = _lbfgsb.types.intvar.dtype
 316
 317     nbd = zeros(n, fortran_int)
 318     low_bnd = zeros(n, float64)
 319     upper_bnd = zeros(n, float64)
 320     bounds_map = {(None, None): 0,
 321                   (1, None): 1,
 322                   (1, 1): 2,
 323                   (None, 1): 3}
 324     for i in range(0, n):
 325         l, u = bounds[i]
 326         if l is not None:
 327             low_bnd[i] = l
 328             l = 1
 329         if u is not None:
 330             upper_bnd[i] = u
 331             u = 1
 332         nbd[i] = bounds_map[l, u]
 333
 334     if not maxls > 0:
 335         raise ValueError('maxls must be positive.')
 336
 337     x = array(x0, float64)
 338     f = array(0.0, float64)
 339     g = zeros((n,), float64)
 340     wa = zeros(2*m*n + 5*n + 11*m*m + 8*m, float64)
 341     iwa = zeros(3*n, fortran_int)
 342     task = zeros(1, 'S60')
 343     csave = zeros(1, 'S60')
 344     lsave = zeros(4, fortran_int)
 345     isave = zeros(44, fortran_int)
 346     dsave = zeros(29, float64)
 347
 348     task[:] = 'START'
 349
 350     n_iterations = 0
 351
 352     while 1:
 353         # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \
 354         _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,
 355                        pgtol, wa, iwa, task, iprint, csave, lsave,
 356                        isave, dsave, maxls)
 357         task_str = task.tobytes()
 358         if task_str.startswith(b'FG'):
 359             # The minimization routine wants f and g at the current x.
 360             # Note that interruptions due to maxfun are postponed
 361             # until the completion of the current minimization iteration.
 362             # Overwrite f and g:
 363             f, g = func_and_grad(x)
 364             if sf.nfev > maxfun:
 365                 task[:] = ('STOP: TOTAL NO. of f AND g EVALUATIONS '
 366                            'EXCEEDS LIMIT')
 367         elif task_str.startswith(b'NEW_X'):
 368             # new iteration
 369             n_iterations += 1
 370             if callback is not None:
 371                 callback(np.copy(x))
 372
 373             if n_iterations >= maxiter:
 374                 task[:] = 'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'
 375             elif sf.nfev > maxfun:
 376                 task[:] = ('STOP: TOTAL NO. of f AND g EVALUATIONS '
 377                            'EXCEEDS LIMIT')
 378         else:
 379             break
 380
 381     task_str = task.tobytes().strip(b'\x00').strip()
 382     if task_str.startswith(b'CONV'):
 383         warnflag = 0
 384     elif sf.nfev > maxfun or n_iterations >= maxiter:
 385         warnflag = 1
 386     else:
 387         warnflag = 2
 388
 389     # These two portions of the workspace are described in the mainlb
 390     # subroutine in lbfgsb.f. See line 363.
 391     s = wa[0: m*n].reshape(m, n)
 392     y = wa[m*n: 2*m*n].reshape(m, n)
 393
 394     # See lbfgsb.f line 160 for this portion of the workspace.
 395     # isave(31) = the total number of BFGS updates prior the current iteration;
 396     n_bfgs_updates = isave[30]
 397
 398     n_corrs = min(n_bfgs_updates, maxcor)
 399     hess_inv = LbfgsInvHessProduct(s[:n_corrs], y[:n_corrs])
 400
 401     task_str = task_str.decode()
 402     return OptimizeResult(fun=f, jac=g, nfev=sf.nfev,
 403                           njev=sf.ngev,
 404                           nit=n_iterations, status=warnflag, message=task_str,
 405                           x=x, success=(warnflag == 0), hess_inv=hess_inv)
 406
 407
 408 class LbfgsInvHessProduct(LinearOperator):
 409     """Linear operator for the L-BFGS approximate inverse Hessian.
 410
 411     This operator computes the product of a vector with the approximate inverse
 412     of the Hessian of the objective function, using the L-BFGS limited
 413     memory approximation to the inverse Hessian, accumulated during the
 414     optimization.
 415
 416     Objects of this class implement the ``scipy.sparse.linalg.LinearOperator``
 417     interface.
 418
 419     Parameters
 420     ----------
 421     sk : array_like, shape=(n_corr, n)
 422         Array of `n_corr` most recent updates to the solution vector.
 423         (See [1]).
 424     yk : array_like, shape=(n_corr, n)
 425         Array of `n_corr` most recent updates to the gradient. (See [1]).
 426
 427     References
 428     ----------
 429     .. [1] Nocedal, Jorge. "Updating quasi-Newton matrices with limited
 430        storage." Mathematics of computation 35.151 (1980): 773-782.
 431
 432     """
 433
 434     def __init__(self, sk, yk):
 435         """Construct the operator."""
 436         if sk.shape != yk.shape or sk.ndim != 2:
 437             raise ValueError('sk and yk must have matching shape, (n_corrs, n)')
 438         n_corrs, n = sk.shape
 439
 440         super().__init__(dtype=np.float64, shape=(n, n))
 441
 442         self.sk = sk
 443         self.yk = yk
 444         self.n_corrs = n_corrs
 445         self.rho = 1 / np.einsum('ij,ij->i', sk, yk)
 446
 447     def _matvec(self, x):
 448         """Efficient matrix-vector multiply with the BFGS matrices.
 449
 450         This calculation is described in Section (4) of [1].
 451
 452         Parameters
 453         ----------
 454         x : ndarray
 455             An array with shape (n,) or (n,1).
 456
 457         Returns
 458         -------
 459         y : ndarray
 460             The matrix-vector product
 461
 462         """
 463         s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho
 464         q = np.array(x, dtype=self.dtype, copy=True)
 465         if q.ndim == 2 and q.shape[1] == 1:
 466             q = q.reshape(-1)
 467
 468         alpha = np.empty(n_corrs)
 469
 470         for i in range(n_corrs-1, -1, -1):
 471             alpha[i] = rho[i] * np.dot(s[i], q)
 472             q = q - alpha[i]*y[i]
 473
 474         r = q
 475         for i in range(n_corrs):
 476             beta = rho[i] * np.dot(y[i], r)
 477             r = r + s[i] * (alpha[i] - beta)
 478
 479         return r
 480
 481     def todense(self):
 482         """Return a dense array representation of this operator.
 483
 484         Returns
 485         -------
 486         arr : ndarray, shape=(n, n)
 487             An array with the same shape and containing
 488             the same data represented by this `LinearOperator`.
 489
 490         """
 491         s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho
 492         I = np.eye(*self.shape, dtype=self.dtype)
 493         Hk = I
 494
 495         for i in range(n_corrs):
 496             A1 = I - s[i][:, np.newaxis] * y[i][np.newaxis, :] * rho[i]
 497             A2 = I - y[i][:, np.newaxis] * s[i][np.newaxis, :] * rho[i]
 498
 499             Hk = np.dot(A1, np.dot(Hk, A2)) + (rho[i] * s[i][:, np.newaxis] *
 500                                                         s[i][np.newaxis, :])
 501         return Hk