python3-bumps 0.7.6-3 / usr / lib / python3 / dist-packages / bumps / quasinewton.py

# Copyright (C) 2009-2010, University of Maryland
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/ or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

# Author: Ismet Sahin
"""
BFGS quasi-newton optimizer.

All modules in this file are implemented from the book
"Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by
J.E. Dennis and Robert B. Schnabel (Only a few minor modifications are done).

The interface is through the :func:`quasinewton` function.  Here is an
example call::

    n = 2
    x0 = [-0.9 0.9]'
    fn = lambda p: (1-p[0])**2 + 100*(p[1]-p[0]**2)**2
    grad = lambda p: array([-2*(1-p[0]) - 400*(p[1]-p[0]**2)*p[0], 200*p[1]])
    Sx = ones(n,1)
    typf = 1                       # todo. see what default value is the best
    macheps = eps
    eta = eps
    maxstep = 100
    gradtol = 1e-6
    steptol = 1e-12                # do not let steptol larger than 1e-9
    itnlimit = 1000
    result = quasinewton(fn, x0, grad, Sx, typf,
                         macheps, eta, maxstep, gradtola, steptol, itnlimit)
    print("status code %d"%result['status'])
    print("x_min=%s, f(x_min)=%g"%(str(result['x']),result['fx']))
    print("iterations, function calls, linesearch function calls",
          result['iterations'],result['evals'],result['linesearch_evals'])
"""
from __future__ import print_function

__all__ = ["quasinewton"]

from numpy import inf, sqrt, isnan, isinf, finfo, diag, zeros, ones
from numpy import array, linalg, inner, outer, dot, amax, maximum

STATUS = {
    1: "Gradient < tolerance",
    2: "Step size < tolerance",
    3: "Invalid point in line search",
    4: "Iterations exceeded",
    5: "Max step taken --- function unbounded?",
    6: "User abort",
    7: "Iterations exceeded in line search",
    8: "Line search step size is too small",
    9: "Singular Hessian",
}


def quasinewton(fn, x0=None, grad=None, Sx=None, typf=1, macheps=None, eta=None,
                maxstep=100, gradtol=1e-6, steptol=1e-12, itnlimit=2000,
                abort_test=None, monitor=lambda **kw: True):
    r"""
    Run a quasinewton optimization on the problem.

    *fn(x)* is the cost function, which takes a point x and returns a scalar fx.

    *x0* is the initial point

    *grad* is the analytic gradient (if available)

    *Sx* is a scale vector indicating the typical values for parameters in
    the fitted result. This is used for a variety of things such as setting
    the step size in the finite difference approximation to the gradient, and
    controlling numerical accuracy in calculating the Hessian matrix.  If for
    example some of your model parameters are in the order of 1e-6, then Sx
    for those parameters should be set to 1e-6. Default: [1, ...]

    *typf* is the typical value for f(x) near the minimum.  This is used along
    with gradtol to check the gradient stopping condition.  Default: 1

    *macheps* is the minimum value that can be added to 1 to produce a number
    not equal to 1.  Default: numpy.finfo(float).eps

    *eta* adapts the numerical gradient calculations to machine precision.
    Default: *macheps*

    *maxstep* is the maximum step size in any gradient step, after normalizing
    by *Sx*. Default: 100

    *gradtol* is a stopping condition for the fit based on the amount of
    improvement expected at the next step.  Default: 1e-6

    *steptol* is a stopping condition for the fit based on the size
    of the step. Default: 1e-12

    *itnlimit* is the maximum number of steps to take before stopping.
    Default: 2000

    *abort_test* is a function which tests whether the user has requested
    abort. Default: None.

    *monitor(x,fx,step)* is called every iteration so that a user interface
    function can monitor the progress of the fit.  Default: lambda \*\*kw: True


    Returns the fit result as a dictionary:

    *status* is a status code indicating why the fit terminated.  Turn the
    status code into a string with *STATUS[result.status]*.  Status values
    vary from 1 to 9, with 1 and 2 indicating convergence and the remaining
    codes indicating some form of premature termination.

    *x* is the minimum point

    *fx* is the value fn(x) at the minimum

    *H* is the approximate Hessian matrix, which is the inverse of the
    covariance matrix

    *L* is the cholesky decomposition of H+D, where D is a small correction
    to force H+D to be positive definite.  To compute parameter uncertainty

    *iterations* is the number of iterations

    *evals* is the number of function evaluations

    *linesearch_evals* is the number of function evaluations for line search
    """
    # print("starting QN")
    # If some input parameters are not specified, define default values for them
    # here. First and second parameters fn and x0 must be defined, others may be
    # passed.  If you want to set a value to a parameter, say to typf, make
    # sure all the parameters before this parameter are specified, in this
    # case fn, x0, grad, and Sx if you want to have default values for grad
    # and Sx, for each enter [].
    # important for also computing fcount (function count)
    n = len(x0)
    if x0 is None:
        x0 = zeros(n)

    if grad is None:
        analgrad = 0
    else:
        analgrad = 1

    if Sx is None:
        Sx = ones(n)
        #Sx = x0 + (x0==0.)
    elif len(Sx) != n:
        raise ValueError("sizes of x0 and Sx must be the same")

    if macheps is None:
        # PAK: use finfo rather than macheps
        macheps = finfo('d').eps

    if eta is None:
        eta = macheps

    fcount = 0                    # total function count
    fcount_ls = 0                # funciton count due to line search

    # If analytic gradient is available then fn will return both function
    # value and analytic gradient.  Otherwise, use finite difference method
    # for estimating the gradient
    if analgrad == 1:
        fc = fn(x0)
        gc = grad(x0)
        fcount = fcount + 1
    else:
        fc = fn(x0)
        gc = fdgrad(n, x0, fc, fn, Sx, eta)
        fcount = fcount + n + 1

    # Check if the initial guess is a local minimizer
    termcode = umstop0(n, x0, fc, gc, Sx, typf, gradtol)
    consecmax = 0

    # Value to return if we fail early
    # Approximately x0 is a critical point
    xf = x0
    ff = fc
    H = L = None
    if termcode == 0:
        H = inithessunfac(n, fc, typf, Sx)

    # STEP 9.
    xc = x0

    # Iterate until convergence in the following loop
    itncount = 0
    while termcode == 0:  # todo. increase itncount
        # print("update",itncount)
        itncount = itncount + 1

        # disp(['Iteration = ' num2str(itncount)])
        # Find Newton step sN
        H, L = modelhess(n, Sx, macheps, H)
        # the vector obtained in the middle
        middle_step_v = linalg.solve(L, -gc)
        sN = linalg.solve(L.transpose(), middle_step_v)   # the last step
        if isnan(sN).any():
            # print("H",H)
            # print("L",L)
            # print("v",middle_step_v)
            # print("Sx",Sx)
            # print("gc",gc)
            termcode = 9
            break

        # Perform line search (Alg.6.3.1). todo. put param order as in the book
        # print("calling linesearch",xc,fc,gc,sN,Sx,H,L,middle_step_v)
        # print("linesearch",xc,fc)
        retcode, xp, fp, maxtaken, fcnt \
            = linesearch(fn, n, xc, fc, gc, sN, Sx, maxstep, steptol)
        fcount += fcnt
        fcount_ls += fcnt
        #plot(xp(1), xp(2), 'g.')

        # Evaluate gradient at new point xp
        if analgrad == 1:
            gp = grad(xp)
        else:
            gp = fdgrad(n, xp, fp, fn, Sx, eta)
            fcount = fcount + n

        # Check stopping criteria (alg.7.2.1)
        consecmax = consecmax + 1 if maxtaken else 0
        termcode = umstop(n, xc, xp, fp, gp, Sx, typf, retcode, gradtol,
                          steptol, itncount, itnlimit, consecmax)

        if abort_test():
            termcode = 6

        # STEP 10.6
        # If termcode is larger than zero, we found a point satisfying one
        # of the termination criteria, return from here.  Otherwise evaluate
        # the next Hessian approximation (Alg. 9.4.1).
        if termcode > 0:
            xf = xp                                        # x final
            ff = fp                                        # f final

        elif not monitor(x=xp, fx=fp, step=itncount):
            termcode = 6

        else:
            H = bfgsunfac(n, xc, xp, gc, gp, macheps, eta, analgrad, H)
            xc = xp
            fc = fp
            gc = gp
        # STOPHERE

    result = dict(status=termcode, x=xf, fx=ff, H=H, L=L,
                  iterations=itncount, evals=fcount, linesearch_evals=fcount_ls)
    #print("result",result, steptol, macheps)
    return result

#------------------------------------------------------------------------------
#@author: Ismet Sahin
# Alg. 9.4.1

# NOTE:
# BFCG Hessian update is performed unless the following two conditions hold
#    (i) y'*s < sqrt(macheps)*norm(s)*norm(y)
#    (ii)


def bfgsunfac(n, xc, xp, gc, gp, macheps, eta, analgrad, H):
    s = xp - xc
    y = gp - gc
    temp1 = inner(y, s)
    # ISMET : I added condition of having temp1 != 0
    if temp1 >= sqrt(macheps) * linalg.norm(s) * linalg.norm(y) and temp1 != 0:
        if analgrad == 1:
            tol = eta
        else:
            tol = sqrt(eta)

        # deal with noise levels in y
        skipupdate = 1
        t = dot(H, s)
        temp_logicals = (abs(y - t) >= tol * maximum(abs(gc), abs(gp)))
        if sum(temp_logicals):
            skipupdate = 0

        # do the BFGS update if skipdate is false
        if skipupdate == 0:
            temp2 = dot(s, t)
            H = H + outer(y, y) / temp1 - outer(t, t) / temp2

    return H


#------------------------------------------------------------------------------
'''
@author: Ismet Sahin
'''


def choldecomp(n, H, maxoffl, macheps):
    minl = (macheps) ** (0.25) * maxoffl

    if maxoffl == 0:
        # H is known to be a positive definite matrix
        maxoffl = sqrt(max(abs(diag(H))))

    minl2 = sqrt(macheps) * maxoffl

    # 3. maxadd is the number (R) specifying the maximum amount added to any
    # diagonal entry of Hessian matrix H
    maxadd = 0

    # 4. form column j of L
    L = zeros((n, n))
    for j in range(1, n + 1):
        L[j - 1, j - 1] = H[j - 1, j - 1] - sum(L[j - 1, 0:j - 1] ** 2)
        minljj = 0
        for i in range(j + 1, n + 1):
            L[i - 1, j - 1] = H[j - 1, i - 1] - \
                sum(L[i - 1, 0:j - 1] * L[j - 1, 0:j - 1])
            minljj = max(abs(L[i - 1, j - 1]), minljj)

        # 4.4
        minljj = max(minljj / maxoffl, minl)

        # 4.5
        if L[j - 1, j - 1] > minljj ** 2:
            # normal Cholesky iteration
            L[j - 1, j - 1] = sqrt(L[j - 1, j - 1])
        else:
            # augment H[j-1,j-1]
            if minljj < minl2:
                minljj = minl2    # occurs only if maxoffl = 0

            maxadd = max(maxadd, minljj ** 2 - L[j - 1, j - 1])
            L[j - 1, j - 1] = minljj

        # 4.6
        L[j:n, j - 1] = L[j:n, j - 1] / L[j - 1, j - 1]

    return L, maxadd

#------------------------------------------------------------------------------
# ALGORITHM 5.6.3

# Ismet Sahin

# function g = fdgrad(n, xc, fc, objfunc, sx, eta)
# g = fdgrad(@obj_function1, 2, [1 -1]', 10, [1 1], eps)

# NOTATION:
#    N : Natural number
#    R : Real number
#    Rn: nx1 real vector
#    Rnxm : nxm real matrix

# INPUTS:
#    n  : the dimension of the gradient vector (N)
#    xc : the current point at which the value of gradient is computed (Rn)
#    fc : function value at xc (R)
#    objfunc : a function handle which is used to compute function values
#    Sx : a n-dim vector, jth entry specifies the typical value of jth param.
# (Rn)
#    eta: equals to 1e-DIGITS where DIGITS is an integer specifying the
# number of reliable digits (R)
# OUTPUT:
#    g : the n-dim finite difference gradient vector (Rn)

# NOTES :
#    hj : is the constant specifying the step size in the direction of jth
# coordinate (R)
#    ej : the unit vector, jth column of the identity matrix (Rn)

# COMMENTS:
#--- FIND STEP SIZE hj
#    1.a : sign(x) does not work for us when x = 0 since this makes the step
# size hj zero which is not allowed. (Step size = 0 => gj = inf.)
#    1.b : evaluation of the step size
#    1.c : a trick to reduce error due to finite precision.  The line xc(j) =
# xc(j) + hj is equivalent to xc = xc + hj * ej where ej is the jth column
# of identity matrix.
#
#--- EVALUATE APPR. GRADIENT
# First evaluate function at xc + hj * ej and then estimate jth entry of
# the gradient.


def fdgrad(n, xc, fc, fn, Sx, eta):

    # create memory for gradient
    g = zeros(n)

    sqrteta = sqrt(eta)
    for j in range(1, n + 1):
        #--- FIND STEP SIZE hj
        if xc[j - 1] >= 0:
            signxcj = 1
        else:
            signxcj = -1                # 1.a

        # 1.b
        hj = sqrteta * max(abs(xc[j - 1]), 1 / Sx[j - 1]) * signxcj

        # 1.c
        tempj = xc[j - 1]
        xc[j - 1] = xc[j - 1] + hj
        hj = xc[j - 1] - tempj

        #--- EVALUATE APPR. GRADIENT
        fj = fn(xc)
        # PAK: hack for infeasible region: point the other way
        if isinf(fj):
            fj = fc + hj
        g[j - 1] = (fj - fc) / hj
        # if isinf(g[j-1]):
        #    print("fc,fj,hj,Sx,xc",fc,fj,hj,Sx[j-1],xc[j-1])

        # now reset the current
        xc[j - 1] = tempj

    #print("gradient", g)
    return g


#------------------------------------------------------------------------------
# @author: Ismet Sahin
# Example call:
# H = inithessunfac(2, f, 1, [1 0.1]')

def inithessunfac(n, f, typf, Sx):
    temp = max(abs(f), typf)
    H = diag(temp * Sx ** 2)
    return H


#------------------------------------------------------------------------------

def linesearch(cost_func, n, xc, fc, g, p, Sx, maxstep, steptol):
    """
ALGORITHM 6.3.1

Ismet Sahin

THE PURPOSE

    is to find a step size which yields the new function value smaller than the
    current function value, i.e. f(xc + alfa*p) <= f(xc) + alfa * lambda * g'p

CONDITIONS

    g'p < 0
    alfa < 0.5

NOTATION:
    N : Natural number
    R : Real number
    Rn: nx1 real vector
    Rnxm : nxm real matrix
    Str: a string

INPUTS
    n : dimensionality (N)
    xc : the current point ( Rn)
    fc : the function value at xc (R)
    obj_func : the function handle to evaluate function values (str like :
       '@costfunction1')
    g : gradient (Rn)
    p : the descent direction (Rn)
    Sx : scale factors (Rn)
    maxstep : maximum step size allowed (R)
    steptol : step tolerance in order to break infinite loop in line search (R)

OUTPUTS
    retcode : boolean indicating a new point xp found (0) or not (1)    (N).
    xp : the new point (Rn)
    fp : function value at xp (R)
    maxtaken : boolean (N)

NOTES:
    alfa : is used to prevent function value reductions which are too small.
       Here we'll use a very small number in order to accept very small
       reductions but not too small.
"""

    maxtaken = 0

    # alfa specifies how much function value reduction is allowable.  The
    # smaller the alfa, the smaller the function value reduction we allow.
    alfa = 1e-4

    # the magnitude of the Newton step
    Newtlen = linalg.norm(Sx * p)

    if Newtlen > maxstep:
        # Newton step is larger than the max acceptable step size (maxstep).
        # Make it equal or smaller than maxstep
        p = p * (maxstep / Newtlen)
        Newtlen = maxstep

    initslope = inner(g, p)

    # "Relative length of p as calculated in the stopping routine"
    # rellength = amax(abs(p) / maximum(abs(xc), Sx))    # this was a bug
    rellength = amax(abs(p) / maximum(abs(xc), 1 / Sx))

    minlambda = steptol / rellength

    lambdaM = 1.0

    # In this loop, we try to find an acceptable next point
    # xp = xc + lambda * p by finding an optimal lambda based on one
    # dimensional quadratic and cubic models
    fcount = 0
    while True:                # 10 starts.
        # next point candidate
        xp = xc + lambdaM * p
        if isnan(xp).any():
            #print("nan xp")
            retcode = 1
            xp, fp = xc, fc
            break
        if fcount > 20:
            #print("too many cycles in linesearch",xp)
            retcode = 2
            xp, fp = xc, fc
            break
        # function value at xp
        fp = cost_func(xp)
        #print("linesearch",fcount,xp,xc,lambdaM,p,fp,fc)
        if isinf(fp):
            fp = 2 * fc  # PAK: infeasible region hack
        fcount = fcount + 1
        if fp <= fc + alfa * lambdaM * initslope:
            # satisfactory xp is found
            retcode = 0
            if lambdaM == 1.0 and Newtlen > 0.99 * maxstep:
                maxtaken = 1
            # return from here
            break
        elif lambdaM < minlambda:
            # step length is too small, so a satisfactory xp cannot be found
            #print("step",lambdaM,minlambda,steptol,rellength)
            retcode = 3
            xp, fp = xc, fc
            break
        else:                            # 10.3c starts
            # reduce lambda by a factor between 0.1 and 0.5
            if lambdaM == 1.0:
                # first backtrack with one dimensional quadratic fit
                lambda_temp = -initslope / (2.0 * (fp - fc - initslope))
                #print("L1",lambda_temp)
            else:
                # perform second and following backtracks with cubic fit
                Mt = array([[1.0/lambdaM**2, -1.0/lambda_prev**2],
                            [-lambda_prev/lambdaM**2, lambdaM/lambda_prev**2]])
                vt = array([[fp - fc - lambdaM * initslope],
                            [fp_prev - fc - lambda_prev * initslope]])
                ab = (1.0 / (lambdaM - lambda_prev)) * dot(Mt, vt)
                # a = ab(1) and b = ab(2)
                disc = ab[1, 0] ** 2 - 3.0 * ab[0, 0] * initslope
                #print("Mt,vt,ab,disc",Mt,vt,ab,disc)
                if ab[0, 0] == 0.0:
                    # cubic model turn out to be a quadratic
                    lambda_temp = -initslope / (2.0 * ab[1, 0])
                    #print("L2",lambda_temp)
                else:
                    # the model is a legitimate cubic
                    lambda_temp = (-ab[1, 0] + sqrt(disc)) / (3.0 * ab[0, 0])
                    #print("L3",lambda_temp)

                if lambda_temp > 0.5 * lambdaM:
                    # larger than half of previous lambda is not allowed.
                    lambda_temp = 0.5 * lambdaM
                    #print("L4",lambda_temp)

            lambda_prev = lambdaM
            fp_prev = fp
            if lambda_temp <= 0.1 * lambdaM:
                # smaller than 1/10 th of previous lambda is not allowed.
                lambdaM = 0.1 * lambdaM
            else:
                lambdaM = lambda_temp

            #print('lambda = ', lambdaM)

    # return xp, fp, retcode
    return retcode, xp, fp, maxtaken, fcount


#------------------------------------------------------------------------------

# @author: Ismet Sahin
# ALGORITHM 1.3.1
def machineeps():
    macheps = 1.0
    while (macheps + 1) != 1:
        macheps = macheps / 2

    macheps = 2 * macheps
    return macheps


#------------------------------------------------------------------------------

def modelhess(n, Sx, macheps, H):
    """
@author: Ismet Sahin.
Thanks to Christopher Meeting for his help in converting this module from
Matlab to Python

ALGORITHM 5.5.1

NOTES:
    Currently we are not implementing steps 1, 14, and 15 (TODO)

This function performs perturbed Cholesky decomposition (CD) as if the input
Hessian matrix is positive definite.  The code for perturbed CD resides in
choldecomp.m file which returns the factored lower triangle matrix L and a
number, maxadd, specifying the largest number added to a diagonal element of
H during the CD decomposition.  This function checks if the decomposition is
completed without adding any positive number to the diagonal elements of H,
i.e. maxadd <= 0.  Otherwise, this function adds the least number to the
diagonals of H which makes it positive definite based on maxadd and other
entries in H.
EXAMPLE CALLS::

         A1 =[2     0    2.4
              0     2     0
              2.4     0     3]

         A2 =[2     0    2.5
               0     2     0
              2.5     0     3]

         A3 =[2     0    10
               0     2     0
              10     0     3]
"""

    # SCALING
    scale_needed = 0                        # ISMET uses this parameter
    if sum(Sx - ones(n)) != 0:
        # scaling is requested by the user
        scale_needed = 1
        Dx = diag(Sx)
        Dx_inv = diag(1.0 / Sx)
        H = dot(Dx_inv, dot(H, Dx_inv))

    # STEP I.
    sqrteps = sqrt(macheps)

    # 2-4.
    H_diag = diag(H)
    maxdiag = max(H_diag)
    mindiag = min(H_diag)

    # 5.
    maxposdiag = max(0, maxdiag)

    # 6. mu is the amount to be added to diagonal of H before the
    # Cholesky decomp. If the minimum diagonal is much much smaller than
    # the maximum diagonal element then adjust mu accordingly otherwise mu = 0.
    if mindiag <= sqrteps * maxposdiag:
        mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag
        maxdiag = maxdiag + mu
    else:
        mu = 0

    # 7. maximum of off-diagonal elements of H
    diag_infinite = diag(inf * ones(n))
    maxoff = (H - diag_infinite).max()

    # 8. if maximum off diagonal element is much much larger than the maximum
    # diagonal element of the Hessian H
    if maxoff * (1 + 2 * sqrteps) > maxdiag:
        mu = mu + (maxoff - maxdiag) + 2 * sqrteps * maxoff
        maxdiag = maxoff * (1 + 2 * sqrteps)

    # 9.
    if maxdiag == 0:            # if H == 0
        mu = 1
        maxdiag = 1

    # 10. mu>0 => need to add mu amount to the diagonal elements: H = H + mu*I
    if mu > 0:
        diag_mu = diag(mu * ones(n))
        H = H + diag_mu

    # 11.
    maxoffl = sqrt(max(maxdiag, maxoff / n))

    # STEP II. Perform perturbed Cholesky decomposition H + D = LL' where D is
    # a diagonal matrix which is implicitly added to H if H is not positive
    # definite. Matrix D has only positive elements. The output variable maxadd
    # indicates the maximum number added to a diagonal entry of the Hesian,
    # i.e. the maximum of D. If maxadd is returned 0, then H was indeed pd
    # and L is the resulting factor.
    # 12.
    L, maxadd = choldecomp(n, H, maxoffl, macheps)

    # STEP III.
    # 13. If maxadd <= 0, we are done H was positive definite.
    if maxadd > 0:
        # H was not positive definite
        # print('WARNING: Hessian is not pd. Max number added to H is ',maxadd)
        maxev = H[0, 0]
        minev = H[0, 0]
        for i in range(1, n + 1):
            offrow = sum(abs(H[0:i - 1, i - 1])) + sum(abs(H[i - 1, i:n]))
            maxev = max(maxev, H[i - 1, i - 1] + offrow)
            minev = min(minev, H[i - 1, i - 1] - offrow)

        sdd = (maxev - minev) * sqrteps - minev
        sdd = max(sdd, 0)
        mu = min(maxadd, sdd)
        H = H + diag(mu * ones(n))
        L, maxadd = choldecomp(n, H, 0, macheps)

    if scale_needed:                # todo. this calculation can be done faster
        H = dot(Dx, dot(H, Dx))
        L = dot(Dx, L)

    return H, L


#------------------------------------------------------------------------------
def umstop(n, xc, xp, f, g, Sx, typf, retcode, gradtol, steptol,
           itncount, itnlimit, consecmax):
    """
#@author: Ismet Sahin

ALGORITHM 7.2.1

Return codes:
Note that return codes are nonnegative integers. When it is not zero, there is
a termination condition which is satisfied.
   0 : None of the termination conditions is satisfied
   1 : Magnitude of scaled grad is less than gradtol; this is the primary
       condition. The new point xp is most likely a local minimizer.  If gradtol
       is too large, then this condition can be satisfied easier and therefore
       xp may not be a local minimizer
   2 : Scaled distance between last two points is less than steptol; xp might be
       a local minimizer.  This condition may also be satisfied if step is
       chosen too large or the algorithm is far from the minimizer and making
       small progress
   3 : The algorithm cannot find a new point giving smaller function value than
       the current point.  The current may be a local minimizer, or analytic
       gradient implementation has some mistakes, or finite difference gradient
       estimation is not accurate, or steptol is too large.
   4 : Maximum number of iterations are completed
   5 : The maximum step length maxstep is taken for last ten consecutive
       iterations.  This may happen if the function is not bounded from below,
       or the function has a finite asymptote in some direction, or maxstep is
       too small.
    """

    termcode = 0
    if retcode == 1:
        termcode = 3
    elif retcode == 2:
        termcode = 7
    elif retcode == 3:
        termcode = 8
    elif retcode > 0:
        raise ValueError("Unknown linesearch return code")
    elif max(abs(g) * maximum(abs(xp), 1 / Sx) / max(abs(f), typf)) <= gradtol:
        # maximum component of scaled gradient is smaller than gradtol.
        # TODO: make sure not to use a too large typf value which leads to the
        # satisfaction of this algorithm easily.
        termcode = 1
    elif max(abs(xp - xc) / maximum(abs(xp), 1 / Sx)) <= steptol:
        # maximum component of scaled step is smaller than steptol
        termcode = 2
    elif itncount >= itnlimit:
        # maximum number of iterations are performed
        termcode = 4
    elif consecmax == 10:
        # not more than 10 steps will be taken consecutively.
        termcode = 5

    return termcode


#------------------------------------------------------------------------------
#@author: Ismet Sahin

# This function checks whether initial conditions are acceptable for
# continuing unconstrained optimization

# f : the function value at x0, i.e. f = f(x0),  (R)
# g : the gradient at x0, (Rn)

# termcode = 0 : x0 is not a critical point of f(x), (Z)
# termcode = 1 : x0 is a critical point of f(x), (Z)

# Note that x0 may be a critical point of the function; in this case, it is
# either a local minimizer or a saddle point of the function.  If the Hessian
# at x0 is positive definite than it is indeed a local minimizer.  Instead of
# checking Hessian, we can also restart the driver program umexample from
# another point which is close to x0.  If x0 is the local minimizer, the
# algorithm will approach it.

def umstop0(n, x0, f, g, Sx, typf, gradtol):
    #consecmax = 0
    if max(abs(g) * maximum(abs(x0), 1./Sx)/max(abs(f), typf)) <= 1e-3*gradtol:
        termcode = 1
    else:
        termcode = 0
    return termcode


#------------------------------------------------------------------------------

def example_call():
    print('***********************************')

    # Rosenbrock function
    fn = lambda p: (1 - p[0])**2 + 100*(p[1] - p[0]**2)**2
    grad = lambda p: array([-2*(1 - p[0]) - 400*(p[1] - p[0]**2)*p[0],
                            200*(p[1] - p[0]**2)])
    x0 = array([2.320894, -0.534223])
    # x0 = array([2.0,1.0])

    result = quasinewton(fn=fn, x0=x0, grad=grad)
    #result = quasinewton(fn=fn, x0=x0)

    print('\n\nInitial point x0 = ', x0, ', f(x0) = ', fn(x0))
    for k in sorted(result.keys()):
        print(k, "=", result[k])


if __name__ == "__main__":
    example_call()