function alpha_star = lsa(x,p,alpha_i)

% The following is an implementaiton of Algorithm 3.5, Nocedal and Wright,
% page 60.

% Note: This code is a simple modified version of a code by Davide Taviani,
% https://gist.github.com/Heliosmaster/1043132

% alpha_1 = a_i
% a0 = a_{i-1}

% Initialization of default parameters
% Here, c1 and c2 are the parameters s recommended by Nocedal and Wright
% (p. 62)

c1 = 1e-4;
c2 = 0.9;
maxit = 100;

alpha_0 = 0;                    % 0th steplength is 0
alpha_max = 10*alpha_i;         % Max serch interval is 10*(initial guess)

[f0,g0] = fun(x);               % Function values at \alpha=0.
df0=dot(p,g0);
 
i=1;

while 1
    fold = fun(x+alpha_0*p);     % Function with previous step-length
    [f,g] = fun(x+alpha_i*p);    % Function with current step-length
    df=dot(g,p);
    
    if( (f > f0+c1*alpha_i*df0) || ((i>1) && (f > fold) ))  % Check for SW1: A sufficent decrease in f or a comparison between f and fold)
        alpha_star = zoom_w(x,p,alpha_0,alpha_i,c1,c2);     % Modify alpha: a suitable value is somewhere in [alpha_0,alpha_1]
        return;
    end
    
    if(abs(df) <= -c2*df0)          % Check for SW2 (Curvature condition)
        alpha_star = alpha_i;       % Current step-length satisfies SW2
        return;
    end
    
    if(df >= 0)
        alpha_star = zoom_w(x,p,alpha_i,alpha_0,c1,c2); % Find a suitable step-length in [alpha_1,alpha_0]
        return;
    end
    
    if(i == maxit)              % Break Clause
        disp('Maximum number of iteration for Line Search reached');
        alpha_star = alpha_i;
        return;
    end
    
    % Update for next loop
    i=i+1;
    alpha_0 = alpha_i;
    
    % Here, I am updating alpha_{i+1} using linear interplation based on
    % rho, rho = 0.8;
    
    rho=0.8;
    alpha_i = rho*alpha_0+(1-rho)*alpha_max;
    
    % alpha_1 = min(alpha_max, alpha_1*3);
      
end

end

% *************************************************************************

function alpha_star = zoom_w(x,p,alpha_lo,alpha_hi,c1,c2)

% The following is an implementaiton of Algorithm 3.6, Nocedal and Wright,
% page 61.

% Note: This code is a simple modified version of a code by Davide Taviani,
% https://gist.github.com/Heliosmaster/1043132

maxit = 20;

[f0,g0] = fun(x);
df0=dot(p,g0);

j = 0;
while 1
    % I use bracketing and bisection to estimate the best value of alpha, 
    % meaning the trial step-length is the middle point of [alpha_lo,alpha_hi]
    
    alpha_j = (alpha_lo+alpha_hi)/2;     
    
    [f,g] = fun(x+alpha_j*p);
    df=dot(p,g);
    
    [f_lo] = fun(x+alpha_lo*p);
    
    if ( (f > f0 + c1*alpha_j*df0) || (f >= f_lo) )   % Test for sufficient decrease or comparison with alpha_lo.
        alpha_hi = alpha_j;                           % Narrow the interval between [alpha_lo,alpha_hi] 
    else
        if abs(df) <= -c2*df0               % Curvature condition (SW2) satisfied
            alpha_star = alpha_j;         
            return;
        end
        if df*(alpha_hi-alpha_lo) >= 0
            alpha_hi = alpha_lo;
        end
        alpha_lo = alpha_j;                 % The interval is now [alpha,alpha_lo]
    end

    if j==maxit
        alpha_star = alpha_j;               % Break Clause
        return;
    end
    
    j = j+1;
end

end