% Class 13 - Optimization
% Coded by NFR on 10.8.2019
% 
% Golden Search (bracketed method, 1d)
f1 = @(x) x^4 + 2*x^3+8*x^2+5*x;
root1 = goldmin(f1,-2,1,1)
% A faster algorithm, built in to Matlab
% also for 1D bracketed problems, use fminbnd
root2 = fminbnd(f1,-2,1)
% fminsearch = open method, which can also be used for
% multiple dimension problems
root3 = fminsearch(f1,-.5)
% example of fminsearch with 2 dimensions
f2 = @(x,y) 2*x.^2+3*y.^2-4*x.*y-y-3*x
% Make contour map of this
xvec = linspace(-10,10,1000);
yvec = linspace(-10,10,1000);
% We need all x-y combinations, use meshgrid
[xmat,ymat] = meshgrid(xvec,yvec);
Zmat = f2(xmat,ymat);
contourf(xvec,yvec,Zmat)
colorbar
% We can see graphical that the point is
% close to (3,2) but let's use fminsearch
% to find it!
% To use fminsearch, we need to cast
% f2 into a 'single' variable format
% X(1) = x
% X(2) = y
f3 = @(X) 2*X(1).^2+3*X(2).^2-4*X(1).*X(2)...
    -X(2)-3*X(1);
root4 = fminsearch(f3,[2 3])

hold on
plot(root4(1),root4(2),'rx','MarkerSize',20)
hold off

% Ladder Example (Problem 7.11 book)
alphavec = (45:135);
np = length(alphavec);
% Preallocate memory for answer vector
L_vec = zeros(np,1);
w1 = 2; w2 = 2; %m

for i = 1:np
    alpha = alphavec(i);
    alpha_rad = alpha/180*pi;
    ladfunc = @(theta) w1/sin(theta)...
        +w2/sin(pi-alpha_rad-theta);
    [theta,Leval] = fminsearch(ladfunc,30/180*pi);
    % Store the answer in preallocated vector
    L_vec(i,1) = Leval;
end
plot(alphavec,L_vec)
xlabel('Angle of hallway (deg)')
ylabel('Maximum Length of Ladder')

% Sometimes you need to find a local min or max
% Need to use a solver that constrains the 
% location of the solution
%
% Example with Himmelblau's function
H1 = @(x,y) (x.^2+y-11).^2+(x+y.^2-7).^2;
xvec = linspace(-5,5,200);
yvec = linspace(-5,5,200);
[xmat,ymat] = meshgrid(xvec,yvec);
Zmat = H1(xmat,ymat);
% 3d plot
surfc(xvec,yvec,Zmat)
% 2d plot
contourf(xvec,yvec,Zmat)
colorbar
xlabel('x')
ylabel('y')
% We want to find minimum in lower left corner
% constrain search to [-5,-1] on both
% x and y
% We need to reform function so it is 
% a single variable, here is a fast way to do this
H2 = @(X) H1(X(1),X(2)); % same thing as above
% but a little faster to write out
root6 = fmincon(H2,[-2,-2],[],[],[],[],[-5,-5],[-1,-1])
hold on
plot(root6(1),root6(2),'rx')