% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------------------------------
% Example 6.1: Powell's Method
%  6.2.3
%  Analytical computation of numbers in Section 6.2.3
%------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% management functions
clear  % clear all variable/information in the workspace - use CAUTION
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  % don't report any warnings like divide by zero etc.

%%%   Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x1 x2 alpha	% symbolic variables
f = 3+(x1-1.5*x2)^2 + (x2-2)^2

% computing the derivative
dfdx1 = diff(f,x1);
dfdx2 = diff(f,x2);
grad =  [dfdx1;dfdx2];

A = [1 0;0 1]
%%%  You will have to enter start vector X0 and seacrh direction
%%% the code will determine analytical solution 
%%%  If alpha is negative - reverse the search direction
%%%  2.14931   -0.43793

X = [0.5 0.5]'

for i = 1:2
gradX = subs(grad,{x1,x2},{X(1),X(2)})
f0 = subs(f,{x1,x2},{X(1),X(2)})

S = -(inv(A)*gradX)
X1 = X + alpha*S;
f1 = subs(f,{x1,x2},{X1(1),X1(2)});

dfa = diff(f1,alpha)
ezplot(f1,[0 5]);
hold on
ezplot(dfa,[0 5]);
grid
ddfa = diff(dfa,alpha);

al = solve(dfa)
astar = double(solve(dfa))

pause
%%% solve for alpha * numerically
anum_5 = fsolve(@falpha,0.5)
anum1_0 = fsolve(@falpha,1)
anum1_5 = fsolve(@falpha,1.5)


% if length(al) == 1
%     aval = double(al)
% end
aval = anum1_5
X1s = subs(X1,alpha,aval)
delx = X1s - X
dfas =  subs(dfa,alpha,aval)
ddfas =  subs(ddfa,alpha,aval)
f1s = subs(f,{x1,x2},{X1s(1),X1s(2)})
gradX1s = subs(grad,{x1,x2},{X1s(1),X1s(2)})

Y = gradX1s - gradX
B = Y*Y'/(Y'*delx)
C=(gradX*gradX')/(gradX'*S)
A = A + B + C
X = X1s
gradX = gradX1s
hold off
close
end
% %S1 = [0 1]'
% S1 = -S1
% % S1 = .0533*[0 1]' + 0.0065*[2.7395 0.8222]'
% X1 = X0 + alpha*S1
% fa = subs(f,{x1,x2},{X1(1),X1(2)})
% dfa = diff(fa,alpha)
% astar = double(solve(dfa))
% 
% for i = 1:length(astar)
%     if (astar(i) > 0) & (astar(i) < 3)
%     fas = subs(fa,astar(i))
%     as = astar(i)
%     X1s = subs(X1,as)
%     end
% end
% 
