% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Modified Newtons Method
%  two dimension version
%  Second Order Method
%  6.4 DFP Method
%  Side constraints are not checked
%  P. Venkataraman - March 2008
%--------------------------------------------------
% Example 6.1:
% f = 3*sin(0.5 + 0.25*x1*x2)*cosx1
%************************************
%                          UpperBound_nVar.m
%						   GoldSection_nVar.m
% and the problem m-file:  Example6_1a.m
%************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% management functions
clear  % clear all variable/information in the workspace - use CAUTION
clear global  %  again use caution - clears global information
clc    % position the cursor at the top of the screen
close   %  closes the figure window
format compact  % avoid skipping a line when writing to the command window
warning off  %#ok<WNOFF> % don't report any warnings like divide by zero etc.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Most values are obtained using symbolic computation
%%%   However optimum stepsize is computed numerically using
%%%   Golden Section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% monitor cpu time
starttime = cputime;
tol = 1.0e-08; lowbound = 0; intvl = 1; ntrials = 20;
xx = [2.5 1.0];  % initial guess for the design
niter = 10;   % number of iterations
functname = 'Example6_1a';
Lc = 'r';
nvar = length(xx);   %% number of design values

for j = 1:nvar
    x(j) = xx(j);
end

if (nvar == 2)
    %*************************
    %  plotting contours -
    %  only for two variables
    %  previous generation code is left in place
    %*************************
    % the plot is centered around initial guess
    % with (+-) delx1, delx2 on either side
    % this can be reset by user

    delx1 = 6;
    delx2 = 5;

    x1 = (xx(1)-delx1):0.1:(xx(1)+delx1);
    x2 = (xx(2)-delx2):0.1:(xx(2)+delx2);
    %     x1 = (0:0.1:4);  % for Example 6_1a
    %     x2 = (-1:0.1:3);
    x1len = length(x1);
    x2len = length(x2);
    %     [X1 X2] = meshgrid(x1,x2);
    %     Vfunctname = strcat('V',functname);
    %     fun = feval(Vfunctname,X1,X2);
    for i = 1:x1len;
        for j = 1:x2len;
            x1x2 =[x1(i) x2(j)];
            fun(j,i) = feval(functname,x1x2);
        end
    end
    colormap(summer);
    [c1,h1] = contourf(x1,x2,fun,[20]);
    colorbar
    %clabel(c1); % remove labelling for clarity
    grid
    xlabel('x_1');
    ylabel('x_2');
    funname = strrep(functname,'_','-');
    title(strcat('Modified Newton Method:',funname));

    % note that contour values are problem dependent
    % the range is problem dependent
    %**************************
    % finished plotting contour
    %***************************

end
clear x1 x2



%%% symbolic calculations
syms x1 x2 f g1 g2 h
f = 3+(x1-1.5*x2)^2 + (x2-2)^2

g1 = diff(f,x1)
g2 = diff(f,x2)
h = jacobian([g1;g2],[x1,x2])

%%% store initial values
x
as(1)=0
xs(1,:)=x;
fs(1)=feval('Example6_1a',x)
for i = 1:niter -1
    % compute Hessian and Gradient values at current x
    hval=subs(h,{x1 x2},{x(1) x(2)})  % hressian
    grad = [subs(g1,{x1 x2},{x(1) x(2)}) subs(g2,{x1 x2},{x(1) x(2)})]'

    fprintf('Iteration: %2i\n',i)
    if grad'*grad <= tol
        fprintf('Exiting - FOC met !!!\n')
        break;
    end
    
    %%% calculate search direction
    s = -inv(hval)*grad;
    sv = double(subs(s,{x1 x2},{x(1) x(2)}))  % make it double
    
    %  call Golden Section for optimum step size
    output = GoldSection_nVar('Example6_1a',tol,x, ...
        sv',lowbound,intvl,ntrials);
    if output(1) ==0
        fprintf('Exiting - Cannot Move - stepsize is zero !!!\n')
        break;
    end
    
    as(i+1) = output(1)
    fs(i+1) = output(2)
    xs(i+1,:)=output(1,3:4);
    x(1,:)=output(1,3:4)
    
        %***********
    % draw lines
    %************

    if (nvar == 2)
        line([xs(i,1) xs(i+1,1)],[xs(i,2) xs(i+1,2)],'LineWidth',2, ...
            'Color',Lc)
        itr = int2str(i);
        x1loc = 0.5*(xs(i,1)+xs(i+1,1));
        x2loc = 0.5*(xs(i,2)+xs(i+1,2));
        %text(x1loc,x2loc,itr);
        % writes iteration number on the line
        if strcmp(Lc,'r')
            Lc = 'k';
        else
            Lc = 'r';
        end

        pause(1)
        %***********************
        % finished drawing lines
        %***********************
    end



end


%%% print total time
totaltime = cputime - starttime;
fprintf('\n\nTotal cpu time (s)= %7.4f \n\n',totaltime)
