function [Top_gazelle_fit,Top_gazelle_pos,Convergence_curve]=GOA2_Laplace(SearchAgents_no,Max_iter,lb,ub,dim,fobj, mu, sigma)

% Inputs:
% SearchAgents_no: number of search agents
% Max_iter: maximum number of iterations
% lb: lower bound of the search space
% ub: upper bound of the search space
% dim: dimension of the problem
% fobj: objective function to be optimized
% mu: mean of the Laplace distribution
% sigma: scale parameter of the Laplace distribution

% Outputs:
% Top_gazelle_fit: fitness value of the best solution found
% Top_gazelle_pos: position of the best solution found
% Convergence_curve: convergence curve of the algorithm

Top_gazelle_pos = zeros(1,dim);
Top_gazelle_fit = inf;
Convergence_curve = zeros(1,Max_iter);
stepsize = zeros(SearchAgents_no,dim);
fitness = inf(SearchAgents_no,1);

% initialize search agents
gazelle = initialization(SearchAgents_no,dim,ub,lb);

% set boundaries
Xmin = repmat(ones(1,dim).lb,SearchAgents_no,1);
Xmax = repmat(ones(1,dim).ub,SearchAgents_no,1);

Iter = 0;
PSRs = 0.34;
S = 0.88;
s = rand();

while Iter<Max_iter
% evaluate top gazelle
for i=1:size(gazelle,1)
Flag4ub = gazelle(i,:) > ub;
Flag4lb = gazelle(i,:) < lb;
gazelle(i,:) = (gazelle(i,:).(~(Flag4ub + Flag4lb))) + ub.Flag4ub + lb.Flag4lb;

fitness(i,1) = fobj(gazelle(i,:));

if fitness(i,1) < Top_gazelle_fit
Top_gazelle_fit = fitness(i,1);
Top_gazelle_pos = gazelle(i,:);
end
end

% keep track of fitness values
if Iter == 0
fit_old = fitness;
Prey_old = gazelle;
end

Inx = (fit_old < fitness);
Indx = repmat(Inx, 1, dim);
gazelle = Indx.Prey_old + ~Indx.gazelle;
fitness = Inx.fit_old + ~Inx.fitness;

fit_old = fitness;
Prey_old = gazelle;

% compute step size
Elite = repmat(Top_gazelle_pos,SearchAgents_no,1);
CF = (1 - Iter/Max_iter)^(2Iter/Max_iter);

RL = levyrnd(SearchAgents_no,dim,mu,sigma); % Levy random number vector
RB = randn(SearchAgents_no,dim); % Brownian random number vector

for i=1:size(gazelle,1)
for j=1:size(gazelle,2)
R = rand();
r = rand();
if mod(Iter,2) == 0
mu_u = -1;
else
mu_u = 1;
end

% exploitation
if r > 0.5
stepsize(i,j) = RB(i,j)(Elite(i,j) - RB(i,j)gazelle(i,j));
gazelle(i,j) = gazelle(i,j) + sRstepsize(i,j);

% exploration
else

if i > size(gazelle,1)/2
stepsize(i,j) = RB(i,j)(RL(i,j)Elite(i,j) - gazelle(i,j));
gazelle(i,j) = Elite(i,j) + Smu_uCFstepsize(i,j);
else
stepsize(i,j) = RL(i,j)(Elite(i,j) - RL(i,j)gazelle(i,j));
gazelle(i,j) = gazelle(i,j) + Smu_uRstepsize(i,j);
end

end
end
end

% update top gazelle
for i=1:size(gazelle,1)

Flag4ub = gazelle(i,:) > ub;
Flag4lb = gazelle(i,:) < lb;
gazelle(i,:) = (gazelle(i,:).(~(Flag4ub+Flag4lb))) + ub.Flag4ub + lb.Flag4lb;

fitness(i,1) = fobj(gazelle(i,:));

if fitness(i,1) < Top_gazelle_fit
Top_gazelle_fit = fitness(i,1);
Top_gazelle_pos = gazelle(i,:);
end
end

% update fitness values
if Iter == 0
fit_old = fitness;
Prey_old = gazelle;
end

Inx = (fit_old < fitness);
Indx = repmat(Inx, 1, dim);
gazelle = Indx.Prey_old + ~Indx.gazelle;
fitness = Inx.fit_old + ~Inx.fitness;
fit_old = fitness;
Prey_old = gazelle;

% apply PSRs
if rand() < PSRs
U = rand(SearchAgents_no,dim) < PSRs;
gazelle = gazelle + CF((Xmin + rand(SearchAgents_no,dim).(Xmax - Xmin)).U);
else
r = rand();
Rs = size(gazelle,1);
stepsize = (PSRs(1-r) + r)(gazelle(randperm(Rs),:) - gazelle(randperm(Rs),:));
gazelle = gazelle + stepsize;
end

Iter = Iter + 1;
Convergence_curve(Iter) = Top_gazelle_fit;
end
end

function Positions = initialization(SearchAgents_no,dim,ub,lb)

Boundary_no = size(ub,2);

if Boundary_no == 1
Positions = rand(SearchAgents_no,dim).(ub-lb) + lb;
end

if Boundary_no > 1
for i = 1:dim
ub_i = ub(i);
lb_i = lb(i);
Positions(:,i) = rand(SearchAgents_no,1).(ub_i-lb_i) + lb_i;
end
end
end

function r = levyrnd(varargin)

if nargin == 1
mu = 0;
sigma = 1;
elseif nargin == 2
mu = varargin{1};
sigma = 1;
elseif nargin == 3
mu = varargin{1};
sigma = varargin{2};
end

r = mu - sigmatan(pi(rand(size(varargin{1}, 1), size(varargin{1}, 2)) - 0.5));
end