can you change my ema code to moema. First, it integrates an enhanced crowding distance function inspired by the Non-dominated Sorting Genetic Algorithm (NSGA), enabling more nuanced solution differentiation. Second, it incorporates a Levy flight approach to significantly augment exploration behavior,here is ema code:

function [FVAL,X,BestFVALIter,Trajectories,fitness_history, position_history] = EMA(NPop,T,lb,ub,dim,FITNESSFCN)  
% Evolutionary Mating Algorithm (EMA)
% EMA attempts to solve problems of the following forms:
%  min F(X)  subject to: lb <= X <= ub
%                  
%  FITNESSFCN   - function handle of the fitness function
%  lb           - lower bounds on X
%  ub           - upper bounds on X
%  T            - number of iterations
%  NPop         - size of the population (class size)  
%  X            - optimal solution determined by EMA
%  FVAL         - value of the fitness function at the minima (X)
%  BestFVALIter - the best fintess function value in each iteration
%  pop          - the population at the end of the specified number of iterations 
% preallocation to store the best objective function of every iteration
% and the objective function value of population
BestFVALIter = NaN(T,1);
obj_dad = NaN(NPop/2,1); % divide the population into 2 groups
obj_mum = NaN(NPop/2,1);
fitness_history=zeros(NPop,T);%  to record  results at each iteration
position_history=zeros(NPop,T,dim);
Trajectories=zeros(dim,T);

% Generation of initial population
pop=initialization(NPop,dim,ub,lb);
pop_dad=pop(1:NPop/2,:);  % Eqn (1) in the paper
pop_mum=pop(NPop/2+1:end,:); % Eqn (2) in the paper

%-----------
%  Evaluation of objective function for initialization
%  Can be vectorized 
for p = 1:NPop/2
    obj_dad(p) = FITNESSFCN(pop_dad(p,:));
    obj_mum(p) = FITNESSFCN(pop_mum(p,:));
end
obj=[obj_dad; obj_mum];
[fmin,I]=min(obj);
best_pop = pop(I,:); % best solution so far
        ub1=ub';
        lb1=lb';

for gen = 1: T    % Start iteration
    k1 = randperm(NPop/2) ;
    k2 = randperm(NPop/2);
    k1x = pop_dad(k1,:)  ;
    k2x = pop_mum(k2,:);
    
    for k=1:NPop/2 % 1 : N/2      
         
    %--------Eqn (3) in the paper is transform into Eqn (4)----%
    %Imates = Isex ratio + [varians(Dad)-varians(Mum)] %ratio
    %for mating (opportunity for sexual selection)
        Imates=1+(var(k1x(k,:))-var(k2x(k,:))); % % Eqn (4) in the paper
        p=randn(1,dim);                 % Eqn (6) in the paper
        r=rand();
        if Imates>=0                    % Upper part of Eqn (5) in the paper
            temp1=p.*k1x(k,:);          % Punnet square 
            temp2=(1-p).*k2x(k,:);
            pop_child=(temp1+temp2);  
            K = rand(1,dim) <0.8;       % Eqn (8) in the paper. Cr is set 0.8 (can be tuned accordingly based on problem)             
            % Eqn (7) in the paper        
            pop_child = best_pop.*(1-K)+pop_child.*K ;  %  Impact of best environment to the offspring            
            % Bounding of the solution
            pop_child = max(min(ub1, pop_child),lb1);
            % evaluate child
            obj_child = FITNESSFCN(pop_child) ; 
            % If the solution improves
            if obj_child<=obj_dad(k1(k))
                pop_dad(k1(k),:)=pop_child;
                obj_dad(k1(k))=obj_child;
            end
            if obj_child<=obj_mum(k2(k))
                pop_mum(k2(k),:)=pop_child;
                obj_mum(k2(k))=obj_child;
            end
            if obj_child<=fmin
                best_pop=pop_child;     % Store the best result if current sol better
                fmin=obj_child;
            end
         end
              
         if Imates<0          % Lower part of Eqn (5) in the paper   
             temp1=p.*k2x(k,:); % Punnet square
             temp2=(1-p).*k1x(k,:);
             pop_child=temp1+temp2;
             K = rand(1,dim) <0.8;  % Eqn (8) in the paper.
             % Eqn (7) in the paper   
             pop_child = best_pop.*(1-K)+pop_child.*K ;% Impact of best environment to the offspring
             % Bounding of the solution
             pop_child = max(min(ub1, pop_child),lb1);
             % evaluate child
             obj_child = FITNESSFCN(pop_child) ;     
              % If the solution improves
             if obj_child<=obj_mum(k2(k))
                 pop_mum(k2(k),:)=pop_child;
                 obj_mum(k2(k))=obj_child;
             end
             if obj_child<=obj_dad(k1(k))
                 pop_dad(k1(k),:)=pop_child;
                 obj_dad(k1(k))=obj_child;
             end
             if obj_child<=fmin
                 best_pop=pop_child;     % Store the best result if current sol better
                 fmin=obj_child;
             end
          end  
          
          %-----Exploration process-----%
          % probablity encountering predator is set to 0.2 (can be tuned depend on problem)%
          if r<0.2
              temp1=best_pop;
              pop_child=rand(1,dim).*temp1;  % Eqn (9) in the paper 
              % Bounding of the solution
              pop_child = max(min(ub1, pop_child),lb1);
              % evaluate child
              obj_child = FITNESSFCN(pop_child) ;
              if obj_child<=fmin
                  best_pop=pop_child; % Store the best result if current sol better
                  fmin=obj_child;
              end               
          end      
     end
     
     % combined updated pool dad and mum 

     pop=[pop_dad;pop_mum];
     obj=[obj_dad;obj_mum];
     [fmin,I]=min(obj);             % Find the best the real best solution
     best_pop = pop(I,:);           % Store the best sol
       
     % record the history at each iter
       fitness_history(:,gen)=obj;
       position_history(:,gen,:)=pop;
       Trajectories(:,gen)=best_pop';
     
     % This is not part of the algorithm but is used to keep track of the
     % best solution determined till the current iteration
     % Savingt the Results    
     BestFVALIter(gen) = fmin;    
     disp(['In iteration #', num2str(gen), ' , target''s objective = ', num2str(BestFVALIter(gen))])
end  % End for iterations

% Extracting the best solution
X = best_pop;
FVAL = BestFVALIter(gen);