% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  DFP Method  for unconstrained minimum
%  First Order Method
%  6.3.3  DFP (Davidon-Fletchet-Powell) method
%  Side constraints are not checked
%  P. Venkataraman - March 2008
%--------------------------------------------------
% A function m-file for the DFP 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 ----------------------------------------------
%  DFP(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
% DFP('Example6_1',[0.5 0.5],20, 1e-08, 0,1 ,20)
%
function ReturnValue = DFP(functname, ...
    dvar0,niter,tol,lowbound,intvl,ntrials)
%%%%%%%%%%%%     global    %%%%%%%%%%%%%
global xfp ffp afp convgfp lenXfp A
%         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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% convergence/stopping criteria
e1 = 1.0e-08; e2 = 1.0e-08; e3 = 1.0e-08;
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('DFP:',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
xfp(1,:) = dvar0;
x = dvar0;
gradNaN = 0;
Lc = 'r';
ffp(1) = feval(functname,x); % value of function at start
afp(1)=0;
grad = (gradfunction(functname,x)); % steepest descent

A = eye(nvar);  % initial metric
% uses MATLAB built in identity matrix function

convgfp(1)=grad*grad';
for i = 1:niter-1
    % determine search direction
%     fprintf('iteration number:  '),disp(i)
    s = (-A*grad')'; % s is used afp a row vector

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

    afp(i+1) = output(1);
    ffp(i+1) = output(2);
    for k = 1:nvar
        xfp(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
    convgfp(i+1)=grad*grad';
    % print convergence value
%     fprintf('gradient length squared:  '),disp(convgfp(i+1));
%     fprintf('objective function value:  '),disp(ffp(i+1));
    %***********
    % draw lines
    %************

    if (nvar == 2)
        line([xfp(i,1) xfp(i+1,1)],[xfp(i,2) xfp(i+1,2)],'LineWidth',2, ...
            'Color',Lc)
        itr = int2str(i);
        x1loc = 0.5*(xfp(i,1)+xfp(i+1,1));
        x2loc = 0.5*(xfp(i,2)+xfp(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 - xfp(i,:);  % difference in the design vector
    lenXfp(i+1) = sqrt(sum(delx.^2));  % length of difference
    if(convgfp(i+1)<= tol)
        fprintf('Kuhn-Tucker Conditions met\n');
        break; end; % convergence criteria
    if(lenXfp(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

    % update the metric here
    % the semicolon hafp been added for Example 6.4
    delx = (x - xfp(i,:))';
    Y = (grad - gradfunction(functname,xfp(i,:)))'; % column vector
    Z = A*Y;
    B = (delx*delx')/(delx'*Y);
    if (norm(Y) <= 1e-08) || (norm(Z) <= 1e-08)
        C = zeros(nvar);
    else C = -(Z*Z')/(Y'*Z);
    end
    A = A + B + C;
end

if gradNaN == 1
    len = length(afp)-1;
else
    len=length(afp);
end
designvar=xfp(len,:);
ReturnValue = [designvar ffp(len)];