%小生境遺傳算法主程序
function f=NicheSGA_1;
%OptimizationStyle='Min';VarietyNumber=1;Xmin=[0];Xmax=[10*pi];PopulationSize=100;GenotypeLength=20;MaxGeneration=100;L=5;
OptimizationStyle='Max';VarietyNumber=1;Xmin=[0];Xmax=[1];PopulationSize=100;GenotypeLength=20;MaxGeneration=100;L=0.18;
Pc=0.9;Pm=0.01;N=PopulationSize/2;Penalty=10^(-30);
initpop=randint(PopulationSize,VarietyNumber*GenotypeLength+VarietyNumber+2);
oldpop=decode(OptimizationStyle,VarietyNumber,PopulationSize,Xmin,Xmax,GenotypeLength,initpop);
[B,index]=sort(oldpop(:,VarietyNumber*GenotypeLength+VarietyNumber+2));
oldpop=[oldpop(index(:),:)];
temp1=oldpop(PopulationSize-N+1:PopulationSize,:);
k=1;
while k<=MaxGeneration
    newpop=selection(VarietyNumber,PopulationSize,GenotypeLength,oldpop);
    randorder=randperm(PopulationSize);
    for i=1:2:(PopulationSize-1)
        parent1=newpop(randorder(i),:);
        parent2=newpop(randorder(i+1),:);
        [child1,child2]=crossover(VarietyNumber,GenotypeLength,parent1,parent2,Pc);
        newpop(randorder(i),:)=child1;
        newpop(randorder(i+1),:)=child2;
    end
    for i=1:PopulationSize
        newpop(i,:)=mutation(VarietyNumber,GenotypeLength,Pm,newpop(i,:));
    end
    newpop=decode(OptimizationStyle,VarietyNumber,PopulationSize,Xmin,Xmax,GenotypeLength,newpop);
    faveragevalue(k)=mean(newpop(:,VarietyNumber*GenotypeLength+VarietyNumber+1));
    if (OptimizationStyle=='Max')
        optimization(k,VarietyNumber+1)=max(newpop(:,VarietyNumber*GenotypeLength+VarietyNumber+1));
    else
        optimization(k,VarietyNumber+1)=min(newpop(:,VarietyNumber*GenotypeLength+VarietyNumber+1));
    end
    [tf,index]=ismember(optimization(k,VarietyNumber+1),newpop(:,VarietyNumber*GenotypeLength+VarietyNumber+1));
    optimization(k,1:VarietyNumber)=newpop(index,VarietyNumber*GenotypeLength+1:VarietyNumber*GenotypeLength+VarietyNumber);
    temp2=[temp1' newpop']';
    temp2=niche(PopulationSize,VarietyNumber,GenotypeLength,N,L,Penalty,temp2);
    [B,index]=sort(temp2(:,VarietyNumber*GenotypeLength+VarietyNumber+2));
    temp2=[temp2(index(:),:)];
    temp1=temp2(PopulationSize+1:PopulationSize+N,:);
    oldpop=temp2(N+1:PopulationSize+N,:);
    k=k+1;
end
if (OptimizationStyle=='Max')
    optimizationf=max(optimization(:,VarietyNumber+1))
else
    optimizationf=min(optimization(:,VarietyNumber+1))
end
[tf,index]=ismember(optimizationf,optimization(:,VarietyNumber+1));
optimizationX=optimization(index,1:VarietyNumber)
x=temp1(:,VarietyNumber*GenotypeLength+1);
y=temp1(:,VarietyNumber*GenotypeLength+VarietyNumber+1);
z=0:0.01:Xmax(1);
plot(x,y,'k.',z,functionf(z),'k-');

%選擇(復制)算子程序
function newpop=selection(VarietyNumber,PopulationSize,GenotypeLength,oldpop);
totalfit=sum(oldpop(:,VarietyNumber*GenotypeLength+VarietyNumber+2));
if (totalfit>1e-15)
    prob=oldpop(:,VarietyNumber*GenotypeLength+VarietyNumber+2)/totalfit;
    prob=cumsum(prob);
    randnum=sort(rand(PopulationSize,1));
    oldorder=1;
    neworder=1;
    while neworder<=PopulationSize
          if (randnum(neworder)<prob(oldorder))
              newpop(neworder,:)=oldpop(oldorder,:);
              neworder=neworder+1;
          else 
              oldorder=oldorder+1;
          end
    end
else
    newpop=oldpop;
end

%交叉算子程序
function [child1,child2]=crossover(VarietyNumber,GenotypeLength,parent1,parent2,Pc);
if(rand<Pc)
    cpoint=round(rand*(VarietyNumber*GenotypeLength-2))+1;
    child1=[parent1(1:cpoint) parent2(cpoint+1:VarietyNumber*GenotypeLength+VarietyNumber+2)];
    child2=[parent2(1:cpoint) parent1(cpoint+1:VarietyNumber*GenotypeLength+VarietyNumber+2)];
else
    child1=parent1;
    child2=parent2;
end

%變異算子程序
function child=mutation(VarietyNumber,GenotypeLength,Pm,parent);
if(rand<Pm)
    mpoint=round(rand*(VarietyNumber*GenotypeLength-1))+1;
    child=parent;
    child(mpoint)=1-parent(mpoint);
else
    child=parent;
end

%小生境子程序
function newpop=niche(PopulationSize,VarietyNumber,GenotypeLength,N,L,Penalty,oldpop);
newpop=oldpop;
for i=1:PopulationSize+N-1
    for j=i+1:PopulationSize+N
        s=0;
        for k=1:VarietyNumber
            s=s+(oldpop(i,VarietyNumber*GenotypeLength+k)-oldpop(j,VarietyNumber*GenotypeLength+k))^2;
        end
        d(i,j)=sqrt(s);
        if (d(i,j)<L)
            if (oldpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)>=oldpop(j,VarietyNumber*GenotypeLength+VarietyNumber+2))
                newpop(j,VarietyNumber*GenotypeLength+VarietyNumber+2)=Penalty;
            else
                newpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)=Penalty;
            end
        end
    end
end

%解碼及求適應度程序
function decodedpop=decode(OptimizationStyle,VarietyNumber,PopulationSize,Xmin,Xmax,GenotypeLength,pop);
decodedpop=pop;
for j=1:PopulationSize
for i=1:VarietyNumber
    decodedpop(j,VarietyNumber*GenotypeLength+i)=(2.^(size(pop(j,(((i-1)*GenotypeLength)+1):(i*GenotypeLength)),2)-1:-1:0)*...
    pop(j,(((i-1)*GenotypeLength)+1):(i*GenotypeLength))')'*(Xmax(i)-Xmin(i))/(2.^GenotypeLength-1)+Xmin(i);         
end
end
for i=1:PopulationSize
    decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+1)=functionf(decodedpop(i,(VarietyNumber*GenotypeLength+1):(VarietyNumber*GenotypeLength+VarietyNumber)));
    if (OptimizationStyle=='Max')
        if (decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+1)>0)
            decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)=1+10*decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+1);
        else
            decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)=1;
        end
    else
        if (decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+1)<0)
            decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)=1-10*decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+1);
        else
            decodedpop(i,VarietyNumber*GenotypeLength+VarietyNumber+2)=1;
        end
    end
end

%計算目標函數值程序
function f=functionf(x);
%f=sin(x).*exp(-0.1*x);
%f=sin(5*3.14159*x).^6;
f=exp(-2*log(2)*((x-0.1)/0.8).^2).*sin(5*3.14159*x).^6;