Why does the vectorized function give wrong result?

I have a function Fun1. You converted it into a vectorized form. I called it Fun2. Then I ran both of them with an algorithm "IGWO". Each time I run them, they give different results. The non-vectorized function Fun1 gives correct result but the vectorized function Fun2 gives wrong result. Why it is so?

Details: Inside the main1.m program, I have a desired vector u on line2. I estimate this desired vector with the help of the algorithm IGWO. For that I call the given functions Fun1 or Fun2 inside the algorithm so that at the end I get the result nearly like u. So I get it with Fun1 but Fun2 gives wrong result though Fun2 is a vectorized form of Fun1. Why it is so?

clear all;clc
u=[1 3 33 53];% Desired Vector 
dim=length(u);
lb=0*ones(1,dim);
ub1=10*ones(1,dim/2);
ub2=90*ones(1,dim/2);
ub=[ub1 ub2];
Runs=2;
one=zeros(Runs,1);
time1=zeros(Runs,1);
temp=zeros(Runs,dim);
two=zeros(Runs,dim);
nn=0;
for n=1:Runs %------------(2)
    nn=nn+1;
    [Alpha_score,Alpha_pos,time]=IGWO(dim,300,1000,lb,ub,@(b)Fun1(b,u))
    % [Alpha_score,Alpha_pos,time]=IGWO(dim,300,1000,lb,ub,@(b)Fun2(b,u));
    
    one(nn)=Alpha_score;
    temp(nn,:)=Alpha_pos;
    time1(nn)=time; 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %-- Swapping
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [~, ix] = sort(u); 
    [~, ix1(ix)] = sort(temp(nn,:)); 
    two(nn,:) = temp(nn,ix1); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save the Workspace data 
%%%%%%%%%%%%%%%%%%%%%%%%%%
%save 2sn0 
function e=Fun1(b,u)
[R,C]=size(b);
P=C/2;
M=2*C;
%%%%%%%%%%%%%%%%%%%%
% Swapping vector b
%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); 
[~, ix1(ix)] = sort(b); 
b = b(ix1);
xo=zeros(1,M);
for k=1:M 
    for i=1:P 
        xo(1,k)=xo(1,k)+1*exp(1i*((k-1)*(-pi/2)*sind(u(P+i))+((k-1)^2*pi/(16*u(i)))*cosd(u(P+i))^2));
    end 
end
xe=zeros(1,M);
for k=1:M
    for i=1:P
        xe(1,k)=xe(1,k)+1*exp(1i*((k-1)*(-pi/2)*sind(b(P+i))+((k-1)^2*pi/(16*b(i)))*cosd(b(P+i))^2));
    end
end
%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%
e=norm(xo-xe).^2/(M);
end
function e=Fun2(b,u)
[R,C]=size(b);
P=C/2;
M=2*C;
k = (1:M).';
i = (1:P);
xo = sum(1*exp(1i*((k-1).*(-pi/2).*sind(b(P+i))+((k-1).^2.*pi./(16*b(i))).*cosd(b(P+i)).^2)),2);
xe = sum(1*exp(1i*((k-1).*(-pi/2).*sind(b(P+i))+((k-1).^2.*pi./(16*b(i))).*cosd(b(P+i)).^2)),2);
%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%
e=norm(xo-xe).^2/(M);
end
%This function is used for L-SHADE bound checking
function vi = boundConstraint (vi, pop, lu)
% if the boundary constraint is violated, set the value to be the middle
% of the previous value and the bound
%
[NP, D] = size(pop); % the population size and the problem's dimension
%% check the lower bound
xl = repmat(lu(1, :), NP, 1);
pos = vi < xl;
vi(pos) = (pop(pos) + xl(pos)) / 2;
%% check the upper bound
xu = repmat(lu(2, :), NP, 1);
pos = vi > xu;
vi(pos) = (pop(pos) + xu(pos)) / 2;
end
%___________________________________________________________________%
% An Improved Grey Wolf Optimizer for Solving Engineering %
% Problems (I-GWO) source codes version 1.0 %
% %
% Developed in MATLAB R2018a %
% %
% Author and programmer: M. H. Nadimi-Shahraki, S. Taghian, S. Mirjalili %
% %
% e-Mail: [e-mail], [e-mail], [e-mail] %
% %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: M. H. Nadimi-Shahraki, S. Taghian, S. Mirjalili %
% An Improved Grey Wolf Optimizer for Solving %
% Engineering Problems , Expert Systems with %
% Applicationsins, in press, %
% DOI: 10.1016/j.eswa.2020.113917 %
%___________________________________________________________________%
%___________________________________________________________________%
% Grey Wold Optimizer (GWO) source codes version 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: [e-mail] %
% [e-mail] %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, S. M. Mirjalili, A. Lewis %
% Grey Wolf Optimizer, Advances in Engineering %
% Software , in press, %
% DOI: 10.1016/j.advengsoft.2013.12.007 %
% %
%___________________________________________________________________%
% Improved Grey Wolf Optimizer (I-GWO)
%function [Alpha_score,Alpha_pos,Convergence_curve]=IGWO(dim,N,Max_iter,lb,ub,fobj)
function [Alpha_score,Alpha_pos,time]=IGWO(dim,N,Max_iter,lb,ub,fobj)
tic;% By Me
lu = [lb .* ones(1, dim); ub .* ones(1, dim)];
% Initialize alpha, beta, and delta positions
Alpha_pos=zeros(1,dim);
Alpha_score=inf; %change this to -inf for maximization problems
Beta_pos=zeros(1,dim);
Beta_score=inf; %change this to -inf for maximization problems
Delta_pos=zeros(1,dim);
Delta_score=inf; %change this to -inf for maximization problems
% Initialize the positions of wolves
Positions=initialization(N,dim,ub,lb);
Positions = boundConstraint (Positions, Positions, lu);
% Calculate objective function for each wolf
for i=1:size(Positions,1)
    Fit(i) = fobj(Positions(i,:));
end
% Personal best fitness and position obtained by each wolf
pBestScore = Fit;
pBest = Positions;
neighbor = zeros(N,N);
Convergence_curve=zeros(1,Max_iter);
iter = 0;% Loop counter
%% Main loop
while iter < Max_iter
    for i=1:size(Positions,1)
        fitness = Fit(i);
        
        % Update Alpha, Beta, and Delta
        if fitness<Alpha_score
            Alpha_score=fitness; % Update alpha
            Alpha_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness<Beta_score
            Beta_score=fitness; % Update beta
            Beta_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score
            Delta_score=fitness; % Update delta
            Delta_pos=Positions(i,:);
        end
    end
    
    %% Calculate the candiadate position Xi-GWO
    a=2-iter*((2)/Max_iter); % a decreases linearly from 2 to 0
    
    % Update the Position of search agents including omegas
    for i=1:size(Positions,1)
        for j=1:size(Positions,2)
            
            r1=rand(); % r1 is a random number in [0,1]
            r2=rand(); % r2 is a random number in [0,1]
            
            A1=2*a*r1-a; % Equation (3.3)
            C1=2*r2; % Equation (3.4)
            
            D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
            
            r1=rand();
            r2=rand();
            
            A2=2*a*r1-a; % Equation (3.3)
            C2=2*r2; % Equation (3.4)
            
            D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2
            
            r1=rand();
            r2=rand();
            
            A3=2*a*r1-a; % Equation (3.3)
            C3=2*r2; % Equation (3.4)
            
            D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3
            
            X_GWO(i,j)=(X1+X2+X3)/3; % Equation (3.7)
            
        end
        X_GWO(i,:) = boundConstraint(X_GWO(i,:), Positions(i,:), lu);
        Fit_GWO(i) = fobj(X_GWO(i,:));
    end
    
    %% Calculate the candiadate position Xi-DLH
    radius = pdist2(Positions, X_GWO, 'euclidean'); % Equation (10)
    dist_Position = squareform(pdist(Positions));
    r1 = randperm(N,N);
    
    for t=1:N
        neighbor(t,:) = (dist_Position(t,:)<=radius(t,t));
        [~,Idx] = find(neighbor(t,:)==1); % Equation (11) 
        random_Idx_neighbor = randi(size(Idx,2),1,dim);
        
        for d=1:dim
            X_DLH(t,d) = Positions(t,d) + rand .*(Positions(Idx(random_Idx_neighbor(d)),d)...
                - Positions(r1(t),d)); % Equation (12)
        end
        X_DLH(t,:) = boundConstraint(X_DLH(t,:), Positions(t,:), lu);
        Fit_DLH(t) = fobj(X_DLH(t,:));
    end
    
    %% Selection 
    tmp = Fit_GWO < Fit_DLH; % Equation (13)
    tmp_rep = repmat(tmp',1,dim);
    
    tmpFit = tmp .* Fit_GWO + (1-tmp) .* Fit_DLH;
    tmpPositions = tmp_rep .* X_GWO + (1-tmp_rep) .* X_DLH;
    
    %% Updating
    tmp = pBestScore <= tmpFit; % Equation (13)
    tmp_rep = repmat(tmp',1,dim);
    
    pBestScore = tmp .* pBestScore + (1-tmp) .* tmpFit;
    pBest = tmp_rep .* pBest + (1-tmp_rep) .* tmpPositions;
    
    Fit = pBestScore;
    Positions = pBest;
    
    %%
    iter = iter+1;
    neighbor = zeros(N,N);
    Convergence_curve(iter) = Alpha_score; 
    time=toc;
end
end
%___________________________________________________________________%
% Grey Wold Optimizer (GWO) source codes version 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: [e-mail] %
% [e-mail] %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, S. M. Mirjalili, A. Lewis %
% Grey Wolf Optimizer, Advances in Engineering %
% Software , in press, %
% DOI: 10.1016/j.advengsoft.2013.12.007 %
% %
%___________________________________________________________________%
% This function initialize the first population of search agents
function Positions=initialization(SearchAgents_no,dim,ub,lb)
Boundary_no= size(ub,2); % numnber of boundaries
% If the boundaries of all variables are equal and user enter a signle
% number for both ub and lb
if Boundary_no==1
    Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end
% If each variable has a different lb and ub
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

[Sadiq]

Thanks dear Walter Roberson for your help. But Torsten has already corrected it that now gives the same result.

Answer 1

[Torsten]

Compare Fun1 and Fun2. In Fun1, xo is computed using u while in Fun2, both xo and xe are computed using b. Further, the sorting part is not done in Fun2 - I don't know exactly what it's good for and whether it makes a difference.

Try this function - it will give the same output as Fun1:

function e=Fun2a(b,u)
  [R,C]=size(b);
  P=C/2;
  M=2*C;
  k = (1:M).';
  i = (1:P);
  %%%%%%%%%%%%%%%%%%%%
  % Swapping vector b
  %%%%%%%%%%%%%%%%%%%%
  [~, ix] = sort(u); 
  [~, ix1(ix)] = sort(b);  
  b = b(ix1);
  xo = sum(1*exp(1i*((k-1).*(-pi/2).*sind(u(P+i))+((k-1).^2.*pi./(16*u(i))).*cosd(u(P+i)).^2)),2);
  xe = sum(1*exp(1i*((k-1).*(-pi/2).*sind(b(P+i))+((k-1).^2.*pi./(16*b(i))).*cosd(b(P+i)).^2)),2);
  %%%%%%%%%%%%%%%%%%
  % MSE
  %%%%%%%%%%%%%%%%%%
  e=norm(xo-xe).^2/(M);
end

[Sadiq]

Thanks a lot dear Torsten for your kind help. Yes, now it gives the same result.