%  本程序是基于人工免疫原理的天氣預報系統的MATLAB實現程序
%  用于免疫系統的學習訓練
%  具體算法請參考該文獻 
%  易世華2006年3月15日
tic                                     %計時開始
clc
clear  
 load data1
 
 a=10;                                  %定義克隆率
 beta=0.2;                              %定義變異概率
 resource=500;                          %定義系統總資源
 sita=0.6;                              %定義刺激閾值
 
 agnor=1/sqrt(12)*ag(:,1:12);
 n=size(ag,1);
 
 %計算整個抗原集合的親和力平均值d
 d=0;  
 for i=1:n-1
     for j=i+1:n
         d=d+sqrt((agnor(i,:)-agnor(j,:))*(agnor(i,:)-agnor(j,:))');
     end
 end
 d=d*2/(n*(n-1))

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %免疫學習系統
 B=zeros(50,17,8);                      %初始化B細胞集合
 M=zeros(15,17,8);                      %初始化M細胞集合
 nb=zeros(1,8);
 nm=nb;
 zero=zeros(1,17);
 
 for i=1:n                              %第i個抗原     
     
     %(1)克隆選擇過程     
     c=ag(i,13)*4+ag(i,14)*2+ag(i,15)+1;%計算抗原(B、M細胞)的類別c
     if  sum(sum(M(:,1:15,c)))
         for j=1:nm(c)
             M(j,16,c)=1-sqrt((M(j,1:12,c)-agnor(i,:))*(M(j,1:12,c)-agnor(i,:))');
         end
         [m num]=max(M(:,16,c));
         mclone=M(num,1:15,c);
         clonen=a*m;
     else
         mclone=[agnor(i,:),ag(i,13:15)];
         nm(c)=nm(c)+1;
         M(nm(c),1:15,c)=mclone;
         clonen=a;
     end
     for j=1:clonen                     %對記憶細胞進行克隆
         rate=rand;
         if rate>beta
            muta=mutate(mclone);        %對抗原的特征向量和類別進行變異
            c2=muta(13)*4+muta(14)*2+muta(15)+1;
            nb(c2)=nb(c2)+1;
            B(nb(c2),1:15,c2)=muta;
            B(nb(c2),16,c2)=1-sqrt((muta(1:12)-agnor(i,:))*(muta(1:12)-agnor(i,:))'); 
            if c2~=c
                B(nb(c2),16,c2)=B(nb(c2),16,c2)*0.8;%與抗原相同類別的B細胞具有更高的刺激值，同時分配更多的資源
            end 
            B(nb(c2),17,c2)=B(nb(c2),16,c2)*a;
        end
     end
     nb(c)=nb(c)+1;
     B(nb(c),1:15,c)=mclone;
     B(nb(c),16,c)=1-sqrt((mclone(1:12)-agnor(i,:))*(mclone(1:12)-agnor(i,:))'); 
     B(nb(c),17,c)=B(nb(c),16,c)*a; 
     
    %(2)動態演化
    si=0;
    while si~=100                        %while    
    for c2=1:8                           %資源競爭算法
        A=sum(B(:,17,c2));
        b1=sum(sum(B(:,16,:)));          %待改進
        b2=resource*sum(B(:,16,c2))/b1;
        while A>b2
            [m,num]=min(B(1:nb(c2),16,c2));            
            B(num,:,c2)=B(nb(c2),:,c2);
            B(nb(c2),:,c2)=zero;
            nb(c2)=nb(c2)-1;
            A=A-m*a;
        end
    end 
    
    for num=1:4                         %隨機選4個B細胞再次進行克隆選擇
        c2=ceil(8*rand);
        while nb(c2)==0
            c2=ceil(8*rand);
        end
        j=ceil(nb(c2)*rand);
        mclone2=B(j,1:15,c2);
        clonen=a*B(j,16,c2);        
        for j=1:clonen
            rate=rand;
            if rate>beta
                muta=mutate(mclone2);
                c2=muta(13)*4+muta(14)*2+muta(15)+1;
                nb(c2)=nb(c2)+1;
                B(nb(c2),1:15,c2)=muta; 
                B(nb(c2),16,c2)=1-sqrt((muta(1:12)-agnor(i,:))*(muta(1:12)-agnor(i,:))');
                if c2~=c
                    B(nb(c2),16,c2)=B(nb(c2),16,c2)*0.8;%與抗原相同類別的B細胞具有更高的刺激值，同時分配更多的資源
                end        
                B(nb(c2),17,c2)=B(nb(c2),16,c2)*a;
            end
        end
    end                
        
    for c2=1:8                          %動態演化過程結束條件                    
        if nb(c2)==0
            continue;
        end
        si=sum(B(:,16,c2))/nb(c2);
        if si>sita             
            si=100;
            break;
        end
    end
    end                                 %while
     
    %(3)生成記憶細胞
    [m,num]=max(B(1:nb(c),16,c));
    mcand=B(num,1:15,c);
    nm(c)=nm(c)+1;
    M(nm(c),1:15,c)=mcand;
    affinity=sqrt((mcand(1:12)-mclone(1:12))*(mcand(1:12)-mclone(1:12))');
    if affinity<d/4
%         [m,num]=min(M(1:nm(c),16,c));            
%         M(num,:,c)=M(nm(c),:,c);
%         M(nm(c),:,c)=zero;
%         nm(c)=nm(c)-1;
        
%         for j=1:nm(c)
%             if sum(M(j,1:15,c)-mclone)==0
%                 M(j,:,c)=B(nm(c),:,c);
%                 M(nm(c),:,c)=zero;
%                 nm(c)=nm(c)-1;
%             end
%         end
    end    
      
  end   
  %免疫學習系統結束
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  save data2 M nm ag agnor n            %保存訓練所得數據
  toc                                   %計時結束