%% minmax2dsolRosenbrock.m  Rosenbrock function
%  Solutions to ex4, part b)
% 
%% Unconstrained Optimization Example
%  "Classical" test function for optimization routines.
%  Banana-shaped valley, min (1,1), value =0, steep walls, shallow bottom
%  near min-point.
% 
% $$f(x,y)=(a-x)^2 + b(y-x^2)^2$$
%
% Take a=1, b=100
% 
% Plot the function to get an idea of where it is minimized

%% 
close all
format compact
a=1;b=100;
f=@(x,y)(a-x).^2 + b*(y-x.^2).^2
fsurf(f,[-1.5,1.5 -1.5 3],'ShowContours','on')
colorbar
% title('f(x,y)=(a-x).^2 + b*(y-x.^2).^2, a=1,b=100')
xlabel('x');ylabel('y')
hold on
plot3(1,1,0,'.r','MarkerSize',15)
%%  Closer to min point
figure
fsurf(f,[0.5,1.5 0.5 1.5],'ShowContours','on')
xlabel('x');ylabel('y')
hold on
plot3(1,1,0,'.r','MarkerSize',15)

%% Contours
figure
clf
x=linspace(-1.5,1.5,200);
y=linspace(-1.5,3,200);
[X,Y]=meshgrid(x,y);
contour(x,y,f(X,Y),60)
title('f(x,y)=(a-x)^2 + b(y-x^2)^2, a=1,b=100')
hold on
contour(x,y,f(X,Y),[0,0.1])
grid on
%% Contours closer to min pt.
figure
x=linspace(0.5,1.5,200);
y=linspace(0.5,1.5,200);
[X,Y]=meshgrid(x,y);
contour(x,y,f(X,Y),60)
hold on
contour(x,y,f(X,Y),[0,0.001 0.005 0.01 0.02 0.1])
grid on
%%
Z=f(X,Y);
min(Z(:))
%% Minimization algorithm: "brute force" vectorization
%  Domain: [-1.5,1.5]x[-1.5,3]  (Previous figure,larger domain)
% 
N=[200 400 800 1600];
Nlen=length(N);
minpts=zeros(3,Nlen);
T=zeros(1,Nlen);
for k=1:Nlen
    tic
    x=linspace(-1.5,1.5,N(k));
    y=linspace(-1.5,3,N(k));
    %x=linspace(0.5,1.5,N(k));y=x; 
    % tic
    [X,Y]=meshgrid(x,y);
    % tic
    Z=f(X,Y);
    % tic
    [minz,minind]=min(Z(:));
    T(k)=toc;
    minpts(:,k)=[X(minind);Y(minind);minz];
end   
format long
minpts
format short
T
ratios=T(2:end)./T(1:end-1)
% 2^nd run (First takes for setup etc., so run at least twice.)
%{
T =
    0.0065    0.0025    0.0089    0.0387
ratios =
    0.3855    3.5368    4.3532
%}


%%
%{
figure
x=linspace(-1.5,1.5,50);
y=linspace(-1.5,3,50);
[X,Y]=meshgrid(x,y);
surf(x,y,f(X,Y))
%}
%% See how fminsearch handles Rosenbrock.
%
% Repeat definitions:
a=1;b=100;
f=@(x,y)(a-x).^2 + b*(y-x.^2).^2
%% Objective function must have intial point in the form of a vector.
%  Just like we had in part a):
%
fun = @(x) f(x(1),x(2));
%%
% One can set optional arguments with |optimset| (remember odeset)
% 
%%
opts=optimset('Display','iter-detailed');
opts.TolFun=1e-10;
opts.TolX=1e-10;
opts.MaxIter=200;
opts.UseParallel=1;  % Presumably doesn't do anything good as
                     % algorithm is purely sequential.
% opts.MaxFunEvals=1000;
%% Call fminsearch with opts set above
tic
[x,fval]=fminsearch(fun,[-1 2],opts);
TfminS=toc
minptfminS=[x fval]
%%
opts.Display
% Iterations: 149
% Accuracy great, because: 
opts.TolFun
opts.TolX
%%
% Timing: Slower than simple vector-min
%