function [MRT,EDGES,cost]=ACA_QoS_MR(C,D,S,E,Dmax,K,M,Alpha,Beta,Gamma,Tau,Rho,Q)
%% Ant Colony Algorithm for QoS Multicast Routing
%  QoS組播路由蟻群算法

%  GreenSim團隊原創作品，轉載請注明

%  此源代碼完整無刪減，請放心使用

%  欲與原作者技術交流請發郵件aihuacheng@gmail.com
%% 輸入參數列表
%  C            費用鄰接矩陣（N×N）
%  D            延時鄰接矩陣（N×N）
%  S            源節點
%  E            組播目的節點（行向量）
%  Dmax         延時約束
%  K            迭代次數（指螞蟻出動多少波）
%  M            螞蟻個數（每一波螞蟻有多少個）
%  Alpha        表征信息素重要程度的參數
%  Beta         表征啟發式因子（費用）重要程度的參數
%  Gamma        表征啟發式因子（延時）重要程度的參數
%  Tau          初始信息素矩陣
%  Rho          信息素蒸發系數
%  Q            信息素增加強度系數
%% 輸出參數列表
%  MRT          最優組播樹（01鄰接矩陣表示）
%  EDGES        最優組播樹所有的邊
%  cost         最優組播樹的費用
%%

%% 第一步：變量初始化
N=size(C,1);%網絡節點個數為N
P=length(E);%目的節點個數為M
MRT=zeros(N,N);
cost=inf;
ROUTES=cell(P,K,M);%用細胞結構存儲到每一個目的節點的每一代的每一只螞蟻的爬行路線
DELAYS=inf*ones(P,K,M);%用三維數組存儲每代每個螞蟻爬行到各個目的節點的延時
COSTS=inf*ones(P,K,M);%用三維數組存儲每代每個螞蟻爬行到各個目的節點的費用
%% 第二步：啟動到P個目的節點的K輪螞蟻覓食活動，每輪派出M只螞蟻
for p=1:P
    Tau=ones(N,N);
    for k=1:K
        for m=1:M
%%        第三步：狀態初始化
            W=S;%當前節點初始化為起始點
            Path=S;%爬行路線初始化
            PD=0;%爬行路線延時初始化
            PC=0;%爬行路線費用初始化
            TABU=ones(1,N);%禁忌表初始化
            TABU(S)=0;%S已經在初始點了，因此要排除
            CC=C;%費用鄰接矩陣備份
            DD=D;%延時鄰接矩陣備份
%%        第四步：下一步可以前往的節點
            DW=DD(W,:);
            DW1=find(DW<inf);
            for j=1:length(DW1)
                if TABU(DW1(j))==0
                    DW(j)=inf;
                end
            end
            LJD=find(DW<inf);%可選節點集
            Len_LJD=length(LJD);%可選節點的個數
%%        覓食停止條件：螞蟻未遇到食物或者陷入死胡同
            while (W~=E(p))&&(Len_LJD>=1)
%%            第五步：轉輪賭法選擇下一步怎么走
                PP=zeros(1,Len_LJD);
                for i=1:Len_LJD
                    PP(i)=(Tau(W,LJD(i))^Alpha)*(C(W,LJD(i))^Beta)*(D(W,LJD(i))^Gamma);
                end
                PP=PP/(sum(PP));%建立概率分布
                Pcum=cumsum(PP);
                Select=find(Pcum>=rand);
                to_visit=LJD(Select(1));%下一步將要前往的節點
%%            第六步：狀態更新和記錄
                Path=[Path,to_visit];%路徑增加
                PD=PD+DD(W,to_visit);%路徑延時累計
                PC=PC+CC(W,to_visit);%路徑費用累計
                W=to_visit;%螞蟻移到下一個節點
                for kk=1:N
                    if TABU(kk)==0
                        CC(W,kk)=inf;
                        CC(kk,W)=inf;
                        DD(W,kk)=inf;
                        DD(kk,W)=inf;
                    end
                end
                TABU(W)=0;%已訪問過的節點從禁忌表中刪除
                DW=DD(W,:);
                DW1=find(DW<inf);
                for j=1:length(DW1)
                    if TABU(DW1(j))==0
                        DW(j)=inf;
                    end
                end
                LJD=find(DW<inf);%可選節點集
                Len_LJD=length(LJD);%可選節點的個數
%%
            end
%%         第七步：記下每一代每一只螞蟻的覓食路線和路線長度
            ROUTES{p,k,m}=Path;
            if Path(end)==E(p)&&PD<=Dmax
                DELAYS(p,k,m)=PD;
                COSTS(p,k,m)=PC;
            else
                DELAYS(p,k,m)=inf;
                COSTS(p,k,m)=inf;
            end
        end
%%     第八步：更新信息素
        Delta_Tau=zeros(N,N);%更新量初始化
        for m=1:M
            if COSTS(p,k,m)<inf&&DELAYS(p,k,m)<Dmax
                ROUT=ROUTES{p,k,m};
                TS=length(ROUT)-1;%跳數
                Cpkm=COSTS(p,k,m);
                for s=1:TS
                    x=ROUT(s);
                    y=ROUT(s+1);
                    Delta_Tau(x,y)=Delta_Tau(x,y)+Q/Cpkm;
                    Delta_Tau(y,x)=Delta_Tau(y,x)+Q/Cpkm;
                end
            end
        end
        Tau=(1-Rho).*Tau+Delta_Tau;%信息素揮發一部分，新增加一部分
    end
end

%% 第九步：整理輸出結果
MINCOSTS=NaN*ones(1,K);
allcost=zeros(1,0);
for k=1:K
    for m=1:M
        COSTkm=COSTS(:,k,m);
        DELAYkm=DELAYS(:,k,m);
        if sum(COSTkm)<inf&&sum(DELAYkm)<inf
            Tree=zeros(N,N);
            for p=1:P
                path=ROUTES{p,k,m};
                RLen=length(path);
                for i=1:(RLen-1)
                    Tree(path(i),path(i+1))=1;
                    Tree(path(i+1),path(i))=1;
                end
            end
            TC=Tree.*C;
            for ii=1:N
                for jj=1:N
                    if isnan(TC(ii,jj))
                        TC(ii,jj)=0;
                    end
                end
            end
            mincost=0.5*sum(sum(TC));
            if mincost<cost
                MINCOSTS(1,k)=mincost;
                MRT=Tree;
                cost=mincost;
            end
            allcost=[allcost,cost];
        end
    end
end
MM=triu(MRT);
T1=find(MM==1);
T2=ceil(T1/N);
T3=mod(T1,N);
EDGES=[T3,T2];
%% 繪收斂曲線
figure(1)
COSTS2=zeros(M,K,P);
DELAYS2=zeros(M,K,P);
for p=1:P
    for k=1:K
        for m=1:M
            if COSTS(p,k,m)<inf
                COSTS2(m,k,p)=COSTS(p,k,m);
                DELAYS2(m,k,p)=DELAYS(p,k,m);
            end
        end
    end
end
LC1=zeros(1,K);
LC2=zeros(1,K);
for k=1:K
    costs=COSTS2(:,k,1);
    delays=DELAYS2(:,k,1);
    pos1=find(costs>0);
    pos2=find(delays>0);
    len1=length(pos1);
    len2=length(pos2);
    LC1(k)=sum(costs)/len1;
    LC2(k)=sum(delays)/len2;
end
plot(LC1,'ko-');
hold on
plot(LC2,'bs-');
legend('費用','延時')
title('路徑的費用延時變化情況')


figure(2)
plot(allcost,'b-')
title('組播樹費用收斂曲線')