% reduced_Newton.m%
% This MATLAB® m-file uses a reduced-Newton algorithm with a
% weak line search to solve a set of non-linear algebraic
% equations.

% (the file has been slightly modified for 10.492 from the version 
% used in course 10.34, available on the MIT OpenCourseWare site.)
%
% The input parameters are :
%
% x0 = a column vector of the initial guess of the unknowns
%
% calc_f = the name of a MATLAB® function that calculates 
% the function vector
%
% calc_Jac = the name of a function that calculates the Jacobian
%
% Options = a data structure containing optional flags
% .max_iter = max # of Newton's method iterations
% .max_iter_LS = max # of weak line search iterations
% .rtol = relative tolerance
% .atol = absolute tolerance
% .step_tol = abs. tolerance below which we switch to full Newton's method
% .verbose = return a trajectory matrices containing the history
%   of the Newton's method iterations
% .use_range = if non-zero, limit the maximum magitude of the full Newton
%   step so that the change in each component is not greater than
%   that in the vector .range
% .range = a vector of the ranges for each of the unknowns. Each component
%   of the Newton step
%
% Param = a data structure containing parameters that are to be passed to
% the calc_f and calc_Jac functions
%
% The output parameters are :
%
% x = the final estimate of the solution
%
% iflag = an integer flag that is 1 for convergence,
% 0 for no convergence, and negative for an error
%
% iter_conv = number of iterations required for convergence
%
% x_traj = a matrix where row # j is the solution estimate at iteration j-1
% f_traj = a matrix where row # j is the function vector at iteration j-1

function [x,iflag,iter_conv,x_traj,f_traj] = ... 
    reduced_Newton(x0,calc_f,calc_Jac,Options,Param); 

% First, signal no convergence. 
iflag = 0;  

% Set number of iterations required for convergence. 
iter_conv = 0; 

% Extract number of state variables. 
Nvar = length(x0); 

% Initialize solution estimate. 
x = x0; 

% Calculate initial function vector. 
f = feval(calc_f,x,Param); 
% f = feval(calc_f,x); 
if(length(f) ~= Nvar) 
iflag = -1; 
error('reduced_Newton: calc_f returns vector of improper length'); 
end

% ensure f is a column vector 
if(size(f,1)~=Nvar)
f = f'; 
end

% Obtain initial norm of the function vector for later convergence tests.
f0_norm_inf = max(abs(f)); 
f_norm_2sq = dot(f,f); 

% Record initial state and function vectors in trajectory. 
count_traj = 1; 
x_traj(count_traj,:) = x'; 
f_traj(count_traj,:) = f'; 

% Set the flag telling us to perform weak line searches. 
i_do_LS = 1; 

% Begin Newton's method iterations 
for iter = 1:Options.max_iter 

    % calculate the Jacobian 
Jac = feval(calc_Jac,x,Param); 
% Jac = feval(calc_Jac,x); 

% Solve the set of linear equations for the full line step 
try 
p = Jac\(-f); 
catch
    iflag = -2;
    error('reduced_Newton: full Newton step calculation error'); 
end

% Now, reduce the magnitude of the Newton step if the user has 
% specified a maximum change allowable for each component.
if(Options.use_range) 

    % Calculate the unit vector lying in the Newton line search 
    % direction. 
p_length = norm(p,2); 
p_unit = p/p_length; 

% Calculate the maximum step in this direction allowable under 
% the condition that each state variable must not change by 
% a magnitude greater than the specified range for that variable. 
step_allow = max(abs(Options.range)); 
for ivar=1:Nvar 
try 
    step_ivar = abs(Options.range(ivar)/p_unit(ivar)); 
    if(step_ivar < step_allow) 
step_allow = step_ivar; 
    end 
end
end
step_allow = min(step_allow,p_length); 
p = p_unit*step_allow;
end

% Begin the weak line search
if(i_do_LS) % perform a weak line search 

    for iter_LS = 0:Options.max_iter_LS 
        iconv_LS = 0; 

        % Calculate fractional step length 
lambda = 2^(-iter_LS); 

% Calculate new solution estimate 
x_new = x + lambda*p; 

% Calculate function at the new solution estimate 
f_new = feval(calc_f,x_new,Param); 
% f_new = feval(calc_f,x_new); 

% Check descent criterion 
f_new_norm_2sq = dot(f_new,f_new); 
if(f_new_norm_2sq <= f_norm_2sq) 
    x = x_new; 
    f = f_new; 
    f_norm_2sq = f_new_norm_2sq; 
    iconv_LS = 1; 
    
break; 
end 
    end

    % If we did not satisfy descent condition, update with final result.
if(~iconv_LS) 
    x = x_new; 
    f = f_new; 
    f_norm_2sq = f_new_norm_2sq; 
end

else % use full Newton step instead

    % Calculate new solution estimate 
x = x + p; 

% Calculate function at the new solution estimate 
f = feval(calc_f,x,Param); 
% f = feval(calc_f,x); 
end

% if in verbose mode, record state and function vectors
if(Options.verbose) 
    count_traj = count_traj + 1; 
    x_traj(count_traj,:) = x'; 
    f_traj(count_traj,:) = f'; 
end

% check for convergence to the solution
f_norm_inf = max(abs(f));
i_conv_rel = 0;
if(f_norm_inf <= Options.rtol*f0_norm_inf)
    i_conv_rel = 1; 
end
i_conv_abs = 0; 
if(f_norm_inf <= Options.atol)
    i_conv_abs = 1; 
end
if((i_conv_rel==1)&(i_conv_abs==1))
    iter_conv = iter; 
    iflag = 1; 
    break; 
end

% Check to see whether need to perform a line search 
% at the next step. 
if(f_norm_inf <= Options.step_tol) 
    i_do_LS = 0; 
else
    i_do_LS = 1; 
end
end 
return; 

% Dr. Barry S. Johnston, Copyright 2004.