% ----------------------------------------------------------------
%            GA Algorithm : by Mohammad Esmaeil Hassanzadeh
%              Project : Rastrigin Function Solving by GA
% ----------------------------------------------------------------
clc
clear all
tic
Xmax=[5.12 5.12];
Xmin=-[5.12 5.12];
n=length(Xmin);
N=280;
Cr=0.1;
Iter_max=100;
miu=0.1;
%------------------- create initial population---- start
% cost=(X-1)^2+(Y-2)^2+(Z-3)^2
for t=1:30
    
    for i=1:1000
        Ipop(i,:)=rand(1,n).*(Xmax-Xmin)+Xmin;
        %         cost(i,1)=sum((Ipop(i,:)-X0).^2);
        cost(i,1)=10*n+sum((Ipop(i,:).^2)-10*cos(2*pi*Ipop(i,:)));
    end
    
    Ipop_mix=[Ipop,cost];   % [X Y Z Cost]
    Ipop_sort=sortrows(Ipop_mix,n+1);  % sort the last column(n+1)
    Best_Solution=Ipop_sort(1,:);
    Ipop_select=Ipop_sort(1:N/2,:); % select N/2 of Population
    %------------------- create initial population---- END
    for iter=1:Iter_max   % global
        Prob=abs(Ipop_select(:,n+1)-Ipop_sort(N/2+1,n+1));
        Prob=Prob/sum(Prob);
        
        Cn(1,1)=Prob(1,1);
        for j=2:N/2
            Cn(j,1)=Prob(j,1)+Cn(j-1,1);
        end
        %--------  choosing fathers & mothers  ----- Start
        for i=1:N/4
            alpha=rand();
            for j=1:N/2
                if Cn(j,1)>=alpha
                    loc_f(i,1)=j;
                    break
                end
            end
        end
        
        for i=1:N/4
            alpha=rand();
            for j=1:N/2
                if Cn(j,1)>=alpha
                    loc_m(i,1)=j;
                    break
                end
            end
        end
        
        Ipop_f=Ipop_select(loc_f,:);
        Ipop_m=Ipop_select(loc_m,:);
        
        %--------  choosing fathers & mothers ------ END
        
        %--------  Mutation  ------ Start
        for i=1:N/4
            beta=rand();
            Xmut_f=beta.*Ipop_f(i,1:n)+(1-beta).*Ipop_m(i,1:n);% bacheye pedar ghaleb
            Xmut_m=beta.*Ipop_m(i,1:n)+(1-beta).*Ipop_f(i,1:n);% bacheye madar ghaleb
            
            for j=1:n
                %-------------Dominant Father
                if Cr>rand()
                    Ipop_chf(i,j)=Ipop_f(i,j);
                else
                    Ipop_chf(i,j)=Xmut_f(1,j);
                end
                
                if Ipop_chf(i,j)<Xmin(j)
                    Ipop_chf(i,j)=Xmin(j)
                end
                
                if Ipop_chf(i,j)>Xmax(j)
                    Ipop_chf(i,j)=Xmax(j)
                end
                %-------------Dominant Mother
                if Cr>rand()
                    Ipop_chm(i,j)=Ipop_m(i,j);
                else
                    Ipop_chm(i,j)=Xmut_m(1,j);
                end
                
                if Ipop_chm(i,j)<Xmin(j)
                    Ipop_chm(i,j)=Xmin(j)
                end
                
                if Ipop_chm(i,j)>Xmax(j)
                    Ipop_chm(i,j)=Xmax(j)
                end
            end
        end
        %----------------- Mutation ---------------- END
        
        for i=1:N/4
            cost_chf(i,1)=10*n+sum((Ipop_chf(i,:).^2)-10*cos(2*pi*Ipop_chf(i,:)));
            cost_chm(i,1)=10*n+sum((Ipop_chm(i,:).^2)-10*cos(2*pi*Ipop_chm(i,:)));
            %
            %             cost_chf(i,1)=sum((Ipop_chf(i,:)-X0).^2);
            %             cost_chm(i,1)=sum((Ipop_chm(i,:)-X0).^2);
        end
        
        clear Ipop_mix
        Ipop_mix=[Ipop_select;Ipop_chf cost_chf;Ipop_chm cost_chm];
        
        %------------  jahesh  ----------------- Start
        Num_mut=max(1,round(N*miu));
        N_mut2=randperm(N)';
        
        Xmax2=(1+miu)*Best_Solution(1,1:n);
        Xmin2=(1-miu)*Best_Solution(1,1:n);
        
        B=[Xmax2;Xmin2];
        Xmin1=min(B);
        Xmax1=max(B);
        % be nazare man niazi be in 2 if nist chon etefagh nemioftad
        for j=1:n
            if Xmax1(1,j)>Xmax(1,j)
                Xmax1(1,j)=Xmax(1,j);
            end
            if Xmin1(1,j)<Xmin(1,j)
                Xmin1(1,j)=Xmin(1,j);
            end
        end
        
        for i=1:Num_mut
            for j=1:15
                Ipop_mut(j,:)=rand(1,n).*(Xmax1-Xmin1)+Xmin1;
                %                 cost_mut(j,1)=sum((Ipop_mut(j,:)-X0).^2);
                cost_mut(j,1)=10*n+sum((Ipop_mut(j,:).^2)-10*cos(2*pi*Ipop_mut(j,:)));
                
            end
            
            Ipop_mix_mut=[Ipop_mut,cost_mut];
            Ipop_mix_sort=sortrows(Ipop_mix_mut,n+1);
            if Ipop_mix(N_mut2(i,1),n+1)>Ipop_mix_sort(1,n+1)
                Ipop_mix(N_mut2(i,1),:)=Ipop_mix_sort(1,:);
            end
        end
        %------------  jahesh  ----------------- END
        
        clear Ipop_sort
        Ipop_sort=sortrows(Ipop_mix,n+1);
        if Best_Solution(1,n+1)>Ipop_sort(1,n+1)
            Best_Solution = Ipop_sort(1,:);
        end
        
        clear Ipop_select
        Ipop_select=Ipop_sort(1:N/2,:);
        
        curve(iter,1)=iter;
        curve(iter,2)=Best_Solution(1,n+1);
        
        BEST(t,:)=[t,Best_Solution(1,n+1),Best_Solution(1,1:n)];
    end
    
end

Best_Solution
% Best
SD=std(BEST(:,2));
fprintf('Standard Deviation : %s  ',SD);
disp('  ')
toc
B2=min(BEST(:,2));
B1=find(BEST(:,2)==B2(1));
Best_answer=BEST(B1(1),2)
W2=max(BEST(:,2));
W1=find(BEST(:,2)==W2(1));
Worst_answer=BEST(W1(1),2)
% -----------------  Plot curves  -----------------
x=(curve(:,1))';
y=(curve(:,2))';
subplot(1,2,1)
plot(x,y,'r','Linewidth',1.5)
grid on
xlabel('Iteration ')
ylabel('cost')
set(gca,'Xcolor',[0.75 0.75 0.75])
set(gca,'Ycolor',[0.75 0.75 0.75])
c_axes = copyobj(gca,gcf);
set(c_axes, 'color', 'none', 'xcolor', 'k', 'xgrid', 'off', 'ycolor','k', 'ygrid','off');
legend('Best Solution')
title('Best Solution')
subplot(1,2,2)
plot(BEST(:,1),BEST(:,2),'b','Linewidth',1.5)
hold on
plot(BEST(:,1),(ones(30,1))*std(BEST(:,2)),'r','Linewidth',1.5)
title('best of 50 times solving')
grid on
xlabel('Iteration ')
ylabel('Best Solution')
set(gca,'Xcolor',[0.75 0.75 0.75])
set(gca,'Ycolor',[0.75 0.75 0.75])
c_axes = copyobj(gca,gcf);
set(c_axes, 'color', 'none', 'xcolor', 'k', 'xgrid', 'off', 'ycolor','k', 'ygrid','off');
legend('Best Solution','Standard Deviation')
set(gcf,'Position',[10 50 1350 640])