%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%   Global Optimization Sketch  -------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

%% Issue with fmin

% Consider the function y = e^{-.1x}*sin(x) defined on x \in [0,\infty).  
% Over its support [0,\infty), this function has a global minimum at
% x=4.613. Additionally, it has many local minima. For example, x=10.896 is
% a local minimizer.

% First, use fmincon to find the global minimum of this function. 

% Initial guess % 
param0 = 4;

% Upper & Lower bounds of search space %
lb = 0;
ub = inf;

fun=@(param)y(param);
[solution,~] = fmincon(fun,param0,lb,ub);
solution 

% The solution is 4.613, so it found the global minimizer. But I sort of
% cheated. I picked an initial guess very close to the global minimizer,
% which I knew beforehand. 

% Let me instead pick a point close to the local minimizer x=10.896
param0 = 10; 
[solution,~] = fmincon(fun,param0,lb,ub);
solution 

% Now, it reports a solution of 10.896. Fmincon, or any of the "fmins,"
% will always report back the first minima it finds, whether it is global
% or local. 

%% Global Minimization Algorithms 

% To avoid local minima, we can use a global optimizer. Simulated annealing
% and particle swarm are two covenient choices. Below is some quick code
% for running both. For both, I'll use initial guess param0=10.896 just to
% show that even with a poor initial guess, these algorithms will still
% find the global minimizer (at least most of the time).
param0 = 10.896; 
fun=@(param)y(param);

% Upper & Lower bounds of search space %
lb = 0;
ub = inf;

%%% Simulated Annealing %%% 
options = optimoptions(@simulannealbnd,'Display','diagnose','InitialTemperature',400);
[solution,~] = simulannealbnd(fun,param0,lb,ub);
solution 

%%% Particle Swarm %%%
x0 = repmat(param0,100,1);
options = optimoptions(@particleswarm,'SwarmSize',100,'Display','iter','InitialSwarmMatrix',x0);
[solution,~] = particleswarm(fun,1,lb,ub);
solution 

% Run this section over and over again and you'll see that most of the
% time, both algorithms find the global minimizer despite me initializing
% it right at a local minimizer. Occasionally it gets stuck in local
% minima. There are a lot of settings that can be tweaked to prevent this
% from happening. 

%% Parallelizing Particle Swarm 

% For computationally-intensive models, particle swarm is preferable to
% simulated annealing. The reason is that PS can be parallelized while SA
% cannot. Below I'll show how to run PS in parallel. 

% # of cpus to be used %
cpu=4; % I set # of cpus = # of function evaluations in each iteration. 
% Note that this number cannot be larger than the # of cpus on your
% computer, unless you're running on longleaf. When
% running on longleaf, it is good to have the # of cpus be a multiple of
% 12.

% Upper & Lower bounds of  search space %
lb = 0;
ub = inf;

% Initial guess % 
% PS doesn't need an initial guess. If the swarm size (# of
% particles on each iteration) is N, then PS will by default use N times 1
% points uniformly spaced over the support of the function to be minimized
% as its initial guess.

% If you would like to manually pick an initial guess, you can first pick a point:
param0 = 7;
% Turn this point into a vector, where each element of the vector will be
% evaluated by a separate particle:
x0=repmat(param0,[cpu,1]); % matrix of initial guesses (one guess per particle)
% Then slightly perturb each guess (so they aren't all evaluating at the
% same point):
x0=x0.*(1+.2*(randi(3,size(x0))-2)); 
% The above line randomly scales up or down (or leaves unchanged) all of
% the elements of x0. 

%parpool('local', cpu); %manually create parallel pool (not usually needed)
fun=@(param)y(param);
warning('off','all');
options = optimoptions(@particleswarm,'UseParallel',true,'SwarmSize',cpu,'Display','iter','InitialSwarmMatrix',x0);
[solution,~] = particleswarm(fun,size(param0,2),lb,ub,options);
solution 

% If you run this over and over again, you'll see that most of the time it
% successfully finds the global minimizer x=4.613. Again, it occasionally 
% settles in local minima. To avoid this, it is good to make the swarm size
% as large as possible. I only used 4 particles in this example, so the coverage 
% of the search space is pretty horrible. In practice you should use as many
% particles as is computationally feasible. 


%% Functions

% y = e^{-.1x}*sin(x)

function obj = y(param)

    obj = exp(-.1*param).*sin(param);
    
end