% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Conjugate Gradient Method  for unconstrained minimum
%  First Order Method
%  6.3.2  Conjugate Gradient (Fletcher-Reeves) method
%  Side constraints are not checked
%  P. Venkataraman - March 2008
%--------------------------------------------------
% A function m-file for the Conjugate Gradient Method
%***************************************************
%***************************************************
% requires:     	       UpperBound_nVar.m
%						   GoldSection_nVar.m
% and the problem m-file:  functname.m
%** gradients are numerically evaluated using first
%** forward difference
%************************************************
% NEW: NEEDS to be TURNED ON BY USER
% DEFAULT :  not used
% For two dimensions it will plot a filled contour
% to take advantage of Array programming
% the method will look for the problem file written
% in array form in  :     Vfunctname.m
%  For Example if functname is "Example6_1.m"
%  The vector form should be in "VExample6_1.m"
%************************************************
%%%%%Usage ----------------------------------------------
%  ConjugateGradient(functname,dvar0,niter,tol,lowbound,intvl,ntrials)
%  functname   :  where specific example can be found
%  dvar0       : initial guess for the design vector
%              : length of this vector defines number of variables
%  niter       : number of random Walk iterations
%  tol         : tolerance for golden section and conjugate gradient
%  lowbound    : for golden section and upper bound calculations
%  intv        :  for upper bound calculations
%  ntrials     : for upper bound calculations
%---------------------------------------------
%%%  Example
% ConjugateGradient('Example6_1',[0.5 0.5],20, 1e-08, 0,1 ,20)
%
function ReturnValue = ConjugateGradient(functname, ...
    dvar0,niter,tol,lowbound,intvl,ntrials)
%%%%%%%%%%%%     global    %%%%%%%%%%%%%
global xcg fcg acg convgcg lenXcg delx1 delx2
%         can be used in calling program to print
%        iterative information if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% management functions
clc    % position the cursor at the top of the screen
clf   %  closes the figure window
format compact  % avoid skipping a line when writing to the command window
warning off  % don't report any warnings like divide by zero etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nvar = length(dvar0); % length of design vector or number of variables
% obtained from start vector
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 = (dvar0(1)-delx1):0.1:(dvar0(1)+delx1);
    x2 = (dvar0(2)-delx2):0.1:(dvar0(2)+delx2);
    %     x1 = (0:0.1:4);  % for Example 6_1
    %     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('Conjugate Gradient:',funname));

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

end

%*********************
%  Numerical Procedure
%*********************
% design vector, alpha , and function value is stored
xcg(1,:) = dvar0;
x = dvar0;
gradNaN = 0;
Lc = 'r';
fcg(1) = feval(functname,x); % value of function at start
acg(1)=0;
s = -(gradfunction(functname,x)); % steepest descent

convgcg(1)=s(1)*s(1)';
for i = 1:niter-1
    % determine search direction
    if (i > 1)
        beta = grad*grad'/(s*s');
        s = -grad + beta*s;
    end

    output = GoldSection_nVar(functname,tol,x, ...
        s,lowbound,intvl,ntrials);

    acg(i+1) = output(1);
    fcg(i+1) = output(2);
    for k = 1:nvar
        xcg(i+1,k)=output(2+k);
        x(k)=output(2+k);
    end
    grad= (gradfunction(functname,x)) ;% steepest descent
    if isnan(grad) == 1
        fprintf('\n***  Gradient does not exist: iteration %2i - EXITING  ***\n',i)
        gradNaN = 1;
        break;
    end
    convgcg(i+1)=grad*grad';
    %***********
    % draw lines
    %************

    if (nvar == 2)
        line([xcg(i,1) xcg(i+1,1)],[xcg(i,2) xcg(i+1,2)],'LineWidth',2, ...
            'Color',Lc)
        itr = int2str(i);
        x1loc = 0.5*(xcg(i,1)+xcg(i+1,1));
        x2loc = 0.5*(xcg(i,2)+xcg(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

    %***************************************
    %  apply stopping criteria
    %****************************************
    delx = x - xcg(i,:);  % difference in the design vector
    lenXcg(i+1) = sqrt(sum(delx.^2));  % length of difference
    if(convgcg(i+1)<= tol)
        fprintf('Kuhn-Tucker Conditions met\n');
        break; end; % convergence criteria
    if(lenXcg(i+1)<= tol)
        fprintf('Exit: Design not changing\n');
        break; end; % convergence criteria
    if(i == niter-1)
        fprintf('Exit: Increase number of Iterations\n');
        break; end; % convergence criteria
end
zoom
if gradNaN == 1
    len = length(acg)-1;

else
    len=length(acg);
    %for kk = 1:nvar

end

%  recover final values to return
designvar=xcg(len,:);
%%% return design information at the last iteration
ReturnValue = [designvar fcg(len)];