function [Destination_fitness,bestPositions,Convergence_curve]=ESMA(N,Max_iter,lb,ub,dim,fobj)
bestPositions=zeros(1,dim);
Destination_fitness=inf;%change this to -inf for maximization problems
AllFitness = inf*ones(N,1);%record the fitness of all slime mold
weight = ones(N,dim);%fitness weight of each slime mold
%Initialize the set of random solutions
X=initialization(N,dim,ub,lb);
Convergence_curve=zeros(1,Max_iter);it=1; %Number of iterations
lb=ones(1,dim).*lb; % lower boundary
ub=ones(1,dim).*ub; % upper boundary
z=0.03; % parameter% Main loop
while it <= Max_iter
for i=1:N
% Check if solutions go outside the search space and bring them back
Flag4ub=X(i,:)>ub;
Flag4lb=X(i,:)<lb;
X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
AllFitness(i) = fobj(X(i,:));
end
%sort the fitness
[SmellOrder,SmellIndex] = sort(AllFitness); %Eq.(3)
bestFitness = SmellOrder(1); %Eq.(5)/ Eq.(10)
worstFitness = SmellOrder(N); %Eq.(6)
bestPositions2=X(SmellIndex(2),:); %Eq.(10)
bestPositions3=X(SmellIndex(3),:); %Eq.(10)
bestPositions4=X(SmellIndex(4),:); %Eq.(10) S=bestFitness-worstFitness+eps; % plus eps to avoid denominator zero
%calculate the fitness weight of each slime mold
for i=1:N
for j=1:dim
if i<=(N/2) %Eq.(4)
weight(SmellIndex(i),j) = 1+rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
else
weight(SmellIndex(i),j) = 1-rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
end
end
end
%update the best fitness value and best position
if bestFitness < Destination_fitness
bestPositions=X(SmellIndex(1),:);
Destination_fitness = bestFitness;
end
avgPositions=(bestPositions+bestPositions2+bestPositions3+bestPositions4)/4; %Eq.(10)
C_pool=[bestPositions; bestPositions2; bestPositions3; bestPositions4; avgPositions]; %Eq.(11)
a = atanh(-(it/Max_iter)+1); %Eq.(8)
b = 1-it/Max_iter; %Eq.(9)
% Update the Position of search agents
for i=1:N
if rand<z %Eq.(12a)
X(i,:) = (ub-lb)*rand+lb;
else
p =tanh(abs(AllFitness(i)-Destination_fitness)); %Eq.(7)
vb = unifrnd(-a,a,1,dim);
vc = unifrnd(-b,b,1,dim);
for j=1:dim
r = rand();
Ceq=C_pool(randi(size(C_pool,1)),:);
A = randi([1,N]); % one positions randomly selected from population
if r<p %Eq.(12b)
X(i,j)=bestPositions(j)+vb(j)*(weight(i,j)*Ceq(j)-X(A,j));
else %Eq.(12c)
X(i,j) = vc(j)*X(i,j);
end
end
end
end
Convergence_curve(it)=Destination_fitness;
it=it+1;
end