%遺傳算法解決tsp問題
function [bestpop,trace]=ga_tsp(d,termops,num,pc,cxops,pm,alpha)
%
%————————————————————————
%[bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%d:距離矩陣
%termops:種群帶數(shù)
%num:每帶染色體的個數(shù)
%pc:交叉概率
%cxops:由于本程序采用單點(diǎn)交叉，交叉點(diǎn)的設(shè)置在本程序中沒有很好的解決，所以本文了采用定點(diǎn)，即第cxops，可以隨機(jī)產(chǎn)生。
%pm:變異概率
%alpha:評價函數(shù)eval(vi)=alpha*(1-alpha).^(i-1).
%bestpop:返回的最優(yōu)種群
%trace:進(jìn)化軌跡
%------------------------------------------------
%####@@@##版權(quán)所有！歡迎廣大網(wǎng)友改正，改進(jìn)！##@@@####
%e-mail:tobysidney33@sohu.com
%####################################################
%
citynum=size(d,2);
n=nargin;
% if n<2
%     disp('缺少變量！！')
%     disp('^_^開個玩笑^_^')
% end
if n<2
    termops=500;
    num=50;
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<3
    num=50;
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<4
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<5
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<6
    pm=0.30;
    alpha=0.10;
end
if n<7
    alpha=0.10;
end
if isempty(cxops)
    cxops=3;
end

[t]=initializega(num,citynum);
for i=1:termops
[l]=f(d,t);
[x,y]=find(l==max(l));
trace(i)=-l(y(1));
bestpop=t(y(1),:);
[t]=select(t,l,alpha);
[g]=grefenstette(t);
[g1]=crossover(g,pc,cxops);
[g]=mutation(g1,pm);  %均勻變異
[t]=congrefenstette(g);
end

%---------------------------------------------------------
function [t]=initializega(num,citynum)
for i=1:num
    t(i,:)=randperm(citynum);
end
%-----------------------------------------------------------
function [l]=f(d,t)
[m,n]=size(t);
for k=1:m
    for i=1:n-1
      l(k,i)=d(t(k,i),t(k,i+1));
    end
      l(k,n)=d(t(k,n),t(k,1));
      l(k)=-sum(l(k,:));
end
%-----------------------------------------------------------
function [t]=select(t,l,alpha)
[m,n]=size(l);
t1=t;
[beforesort,aftersort1]=sort(l,2);%fsort from l to u
for i=1:n
    aftersort(i)=aftersort1(n+1-i);      %change 
end
for k=1:n;
    t(k,:)=t1(aftersort(k),:);
    l1(k)=l(aftersort(k));
end
t1=t;
l=l1;
for i=1:size(aftersort,2)
    evalv(i)=alpha*(1-alpha).^(i-1);
end
m=size(t,1);
q=cumsum(evalv);
qmax=max(q);
for k=1:m
    r=qmax*rand(1);
    for j=1:m
        if j==1&r<=q(1)
            t(k,:)=t1(1,:);
        elseif j~=1&r>q(j-1)&r<=q(j)
            t(k,:)=t1(j,:);
        end
    end
end
%--------------------------------------------------
function [g]=grefenstette(t)
[m,n]=size(t);
for k=1:m
    t0=1:n;
   for i=1:n
       for j=1:length(t0)
           if t(k,i)==t0(j)
              g(k,i)=j;
              t0(j)=[];
              break
           end
       end
    end
end
%-------------------------------------------
function [g]=crossover(g,pc,cxops)
[m,n]=size(g);
ran=rand(1,m);
r=cxops;
[x,ru]=find(ran<pc);
if ru>=2
    for k=1:2:length(ru)-1
       g1(ru(k),:)=[g(ru(k),[1:r]),g(ru(k+1),[(r+1):n])];
       g(ru(k+1),:)=[g(ru(k+1),[1:r]),g(ru(k),[(r+1):n])];
       g(ru(k),:)=g1(ru(k),:);
    end
end
%--------------------------------------------
function [g]=mutation(g,pm)    %均勻變異
[m,n]=size(g);
ran=rand(1,m);
r=rand(1,3);      %dai gai jin
rr=floor(n*rand(1,3)+1);
[x,mu]=find(ran<pm);
for k=1:length(mu)
    for i=1:length(r)
        umax(i)=n+1-rr(i);
        umin(i)=1;
        g(mu(k),rr(i))=umin(i)+floor((umax(i)-umin(i))*r(i));
    end
end
%---------------------------------------------------
function [t]=congrefenstette(g)
[m,n]=size(g);
for k=1:m
    t0=1:n;
   for i=1:n
      t(k,i)=t0(g(k,i));
      t0(g(k,i))=[];
   end
end