%  本程序是基于人工免疫原理的天氣預(yù)報(bào)系統(tǒng)的MATLAB實(shí)現(xiàn)程序
%  用于免疫系統(tǒng)的學(xué)習(xí)訓(xùn)練
%  具體算法請(qǐng)參考該文獻(xiàn) 
%  易世華2006年3月15日
tic                                     %計(jì)時(shí)開(kāi)始
clc
clear  
 load data1
 
 a=10;                                  %定義克隆率
 beta=0.2;                              %定義變異概率
 resource=500;                          %定義系統(tǒng)總資源
 sita=0.6;                              %定義刺激閾值
 
 agnor=1/sqrt(12)*ag(:,1:12);
 n=size(ag,1);
 
 %計(jì)算整個(gè)抗原集合的親和力平均值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))

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %免疫學(xué)習(xí)系統(tǒng)
 B=zeros(50,17,8);                      %初始化B細(xì)胞集合
 M=zeros(15,17,8);                      %初始化M細(xì)胞集合
 nb=zeros(1,8);
 nm=nb;
 zero=zeros(1,17);
 
 for i=1:n                              %第i個(gè)抗原     
     
     %(1)克隆選擇過(guò)程     
     c=ag(i,13)*4+ag(i,14)*2+ag(i,15)+1;%計(jì)算抗原(B、M細(xì)胞)的類別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;
         M(nm(c),16,c)=1;
         M(nm(c),17,c)=a;
         clonen=a;
     end
     for j=1:clonen                     %對(duì)記憶細(xì)胞進(jìn)行克隆
         rate=rand;
         if rate>beta
            muta=mutate(mclone);        %對(duì)抗原的特征向量和類別進(jìn)行變異
            c2=muta(13)*4+muta(14)*2+muta(15)+1;
            nb(c2)=nb(c2)+1;
            B(nb(c2),1:15,c2)=muta;            
        end
     end
     nb(c)=nb(c)+1;
     B(nb(c),1:15,c)=mclone;
          
    %(2)動(dòng)態(tài)演化
    si=0;
    while si~=100                        %while 
        
    for c2=1:8
        for j=1:nb(c2)
            B(j,16,c2)=1-sqrt((B(j,1:12,c2)-agnor(i,:))*(B(j,1:12,c2)-agnor(i,:))'); 
            if c2~=c
                B(j,16,c2)=B(j,16,c2)*0.8;%與抗原相同類別的B細(xì)胞具有更高的刺激值，同時(shí)分配更多的資源
            end 
            B(j,17,c2)=B(j,16,c2)*a;
        end
    end        
        
    for c2=1:8                           %資源競(jìng)爭(zhēng)算法
        A=sum(B(:,17,c2));
        b1=sum(sum(B(:,16,:)));          %待改進(jìn)
        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                         %隨機(jī)選4個(gè)B細(xì)胞再次進(jìn)行克隆選擇
        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細(xì)胞具有更高的刺激值，同時(shí)分配更多的資源
                end        
                B(nb(c2),17,c2)=B(nb(c2),16,c2)*a;
            end
        end
    end                
        
    for c2=1:8                          %動(dòng)態(tài)演化過(guò)程結(jié)束條件                    
        if nb(c2)==0
            continue;
        end
        si=sum(B(:,16,c2))/nb(c2);
        if si>sita             
            si=100;
            break;
        end
    end
    end                                 %while
     
    %(3)生成記憶細(xì)胞
    [m,num]=max(B(1:nb(c),16,c));
    mcand=B(num,1:15,c);
    nm(c)=nm(c)+1;
    M(nm(c),:,c)=B(num,:,c);
    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

%           M(nm(c)-1,:,c)=M(nm(c),:,c);
%           M(nm(c),:,c)=zero;
%           nm(c)=nm(c)-1;          
    end
    
      
  end                                   %for i=1:n 
  
%   for c2=1:8
%       if nm(c2)==0
%           [m,num]=max(B(1:nb(c2),16,c2));
%           nm(c2)=nm(c2)+1;
%           M(nm(c2),:,c2)=B(num,:,c2);
%       end
%   end
  %免疫學(xué)習(xí)系統(tǒng)結(jié)束
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  save data2 M nm ag agnor n            %保存訓(xùn)練所得數(shù)據(jù)
  toc                                   %計(jì)時(shí)結(jié)束