% ga_pks.m
% 本程式根據 張智星教授 以及 王進德 教授程式，以及晚學自己一點小意見，改寫而成。
% 南台電機 趙春棠 Chun-Tang Chao, 2007. 5.27.  No Copyright.
% 利用基因演算法，找 Matlab 中著名的 peaks function 的極值。
% 在 Matlab 中，首先執行 >> peaks  即可看到該函數。當然，本程式可以輕易擴充求任意函數的極值。

clear all; clc;
close all;
M=80;           % 偶數為宜
G=100;
len=10;   
umax= 3;
umin= -3;
% 第1代的所有基因
popu =round(rand(M,2*len))

for k=1:1:G
    time(k)=k;
    % 母代及函數值計算
    for j=1:1:M
        % 一次一個基因
        m=popu(j,:) ;       
        % 基因中的第一個變數(二進位)
        % 轉為介於umin與umax間的實數
        m1=m(1:1:len);
        s=0;    % 不宜用 sum 命名喔
        for i=1:1:len
            s=s+m1(i)*2^(len-i);
        end
        x1(j)=(umax-umin)*s/(2^len-1)+umin;
        % 基因中的第二個變數(二進位)
        % 轉為介於umin與umax間的實數
        m2=m(len+1:1:2*len);
        s=0;
        for i=1:1:len
            s=s+m2(i)*2^(len-i);
        end
        x2(j)=(umax-umin)*s/(2^len-1)+umin;
        % 計算目標函數f(x1,x2)值
        F(j)=peaksfun(x1(j),x2(j));
        % F(j)=100*(x1(j)^2-x2(j))^2+(1-x1(j))^2  ; 

    end
    % 紀錄最佳的結果
    % Elitism: find the best two and keep them.
   [Fs_b, indexF]=sort(F);
    best_set(k,:)=[x1(indexF(M)) x2(indexF(M)) Fs_b(M)];
    best_2= popu([indexF(M) indexF(M-1)],:);
    new_popu([1 2],:)=best_2;
    
    % Fitness 計算
    F_positive=F-min(F);        % fitness, keep it positive
    fitness=F_positive/sum(F_positive);  % Scaled fitness序列，總合為1 喔！
    cum_prob=cumsum(fitness);
    
    % 基因Selection & Crossover 
    % 這是在所有的基因中挑選，所以前面保留的最佳兩組，還是可能被挑選到。
    % 但前面保留的最佳兩組，是不會因為基因交配，而被更動喔。
    for i=2:M/2
       tmp = find(cum_prob - rand > 0);
       parent1 = popu(tmp(1), :);
       parent1_hl=popu(tmp(1), 1:len); parent1_hr = popu(tmp(1), len+1:end); 
       tmp = find(cum_prob - rand > 0);
       parent2 = popu(tmp(1), :);
       parent2_hl=popu(tmp(1), 1:len); parent2_hr = popu(tmp(1), len+1:end);
       
       xp = ceil(rand*(len-1)); % 決定交配點 crossover_point
       % 進行交配，並收錄 交配結果
       new_popu(i*2-1, :) = ... 
            [parent1_hl(1:xp) parent2_hl(xp+1:len) parent1_hr(1:xp) parent2_hr(xp+1:len)];
       new_popu(i*2, :) = ...
            [parent2_hl(1:xp) parent1_hl(xp+1:len) parent2_hr(1:xp) parent1_hr(xp+1:len)];
    end
    
    %== MUTATION (elites are not subject to this.)
    mut_rate=0.1;
    mask = rand(M, len*2) < mut_rate; 
    new_popu=xor(new_popu,mask);
    
    % 還原最佳的基因 (最佳基因不可以被 突變 破壞喔！)
    new_popu([1 2],:)= best_2;
    popu=new_popu;
end

best_set
plot(time,best_set(:,3));
xlabel('times'); ylabel('Best Function Value');
best_set(k,1), best_set(k,2)   % x1, x2 value