%% Find Global or Multiple Local Minima
% Lecture4/myGlobalorMultipleLocalMinima.m
%
%%
% <https://se.mathworks.com/help/gads/example-finding-global-or-multiple-local-minima.html> 
%
%% Function to Optimize
%
% This example illustrates how GlobalSearch finds a global minimum efficiently, 
% and how MultiStart finds many more local minima.
%
% The global minimum is at r = 0, with objective function 0. The function g(r) 
% grows approximately linearly in r, with a repeating sawtooth shape. 
% The function h(t) has two local minima, one of which is global.
%%
format compact
close all
%%
type sawtoothxy
%%
t=linspace(0,2*pi);r=linspace(0,20);
h = cos(2*t - 1/2)/2 + cos(t) + 2;
g = (sin(r) - sin(2*r)/2 + sin(3*r)/3 - sin(4*r)/4 + 4) ...
    .*r.^2./(r+1);
subplot(1,2,1)
plot(t,h)
xlabel('\theta');ylabel('h(\theta)')
subplot(1,2,2)
plot(r,g)
xlabel('r');ylabel('g(r)')
%%
figure
x=linspace(-20,20,30); y=x;
[X,Y]=meshgrid(x,y);
Z = sawtoothxy(X,Y);
surf(X,Y,Z);shg
%% The same with fsurf
figure
fsurf(@(x,y) sawtoothxy(x,y), [-20 20])
xlabel('x');ylabel('y')
%% Create the problem structure.
%  Use the 'sqp'-algorithm.
problem = createOptimProblem('fmincon',...
    'objective',@(x)sawtoothxy(x(1),x(2)),...
    'x0',[100,-50],'options',...
    optimoptions(@fmincon,'Algorithm','sqp','Display','off'));
problem
% The start point is [100,-50] 
%% 3. Validate the problem structure by running fmincon:
%
[x,fval] = fmincon(problem)
%% 4. Create the GlobalSearch object, and set iterative display:
gs = GlobalSearch('Display','iter')
%% Run the solver
rng(14,'twister') % for reproducibility
[x,fval] = run(gs,problem)
% The solver found three local minima, and it found the global minimum ...
% near [0,0]
%% Multiple Local Minima Via MultiStart
% (1) sawtoothxy, as above.
% (2) Create the problem structure. Use the fminunc solver with the 
% Algorithm option set to 'quasi-newton'.The reasons for these choices are:
% - unconstrained. (see desicion table)
% - fminunc requires a grad => 'quasi-Newton'
problem = createOptimProblem('fminunc',...
    'objective',@(x)sawtoothxy(x(1),x(2)),...
    'x0',[100,-50],'options',...
    optimoptions(@fminunc,'Algorithm','quasi-newton','Display','off'))
% - Validate the problem structure by running it:
%%
[x,fval] = fminunc(problem)
%% 4. Create a default MultiStart object
ms = MultiStart 
%% 5. Run the solver for 50 iterations, recording the local minima:
[x,fval,eflag,output,manymins] = run(ms,problem,50)

%MultiStart started fminunc from start points with components uniformly 
% distributed between –1000 and 1000. fminunc often got stuck in one of 
% the many local minima. fminunc exceeded its iteration limit or function 
% evaluation limit 40 times.
%% Parallel MultiStart
%
% <https://se.mathworks.com/help/gads/example-parallel-multistart.html>
%
type sawtoothxy
% 
problem = createOptimProblem('fminunc',...
    'objective',@(x)sawtoothxy(x(1),x(2)),...
    'x0',[100,-50],'options',...
    optimoptions(@fminunc,'Algorithm','quasi-newton'));
%%
[x,fval] = fminunc(problem)
%%
ms = MultiStart('UseParallel',true,'Display','iter');
%%
parpool
%%
[x,fval,eflag,output,manymins] = run(ms,problem,50);
%%
Xloc=[manymins.X]
Floc=[manymins.Fval]
size(Xloc)
size(Floc)



