%  本程序是基于人工免疫原理的天氣預(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);                          %找出ag的行數(shù)
 
 %計(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,:))');%可用norm(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ì)胞)的類(lèi)別c
     if  sum(sum(M(:,1:15,c)))          %可用nm(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ì)抗原的特征向量和類(lèi)別進(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;%與抗原相同類(lèi)別的B細(xì)胞具有更高的刺激值，同時(shí)分配更多的資源
            end 
            B(j,17,c2)=B(j,16,c2)*a;    %是否可考慮將原先分配的資源數(shù)疊加進(jìn)去 ？B(j,17,c2)=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;                    % 可改為A=A-B(num,17,c2); 放至[m,num]=min(B(1:nb(c2),16,c2)); 之后
        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;%與抗原相同類(lèi)別的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));   %刪除刺激值最小的記憶細(xì)胞       
%         M(num,:,c)=M(nm(c),:,c);
%         M(nm(c),:,c)=zero;
%         nm(c)=nm(c)-1;
        
        for j=1:nm(c)                   %刪除記憶細(xì)胞mclone,識(shí)別效果   
            if sum(M(j,1:15,c)==mclone)==15
                M(j,:,c)=M(nm(c),:,c);
                M(nm(c),:,c)=zero;
                nm(c)=nm(c)-1;
                break
            end
        end

%           M(nm(c)-1,:,c)=M(nm(c),:,c); %刪除前一個(gè)記憶細(xì)胞
%           M(nm(c),:,c)=zero;
%           nm(c)=nm(c)-1;          
    end
    
      
  end                                   %for i=1:n 
  
%   for c2=1:8        % 如果某類(lèi)記憶細(xì)胞個(gè)數(shù)為零，則加進(jìn)一個(gè)B細(xì)胞，
%       if nm(c2)==0  % 這樣做沒(méi)多少道理，因?yàn)槟愁?lèi)無(wú)訓(xùn)練樣本，則顯然不會(huì)有記憶細(xì)胞
%           [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é)束