function [Ab] = optd1()%(ts,N,Nc,beta,gen)

% function [M] = optd1(ts,N,Nc,beta,gen)
% M     -> matrix of memory cells  記憶細胞矩陣
% ts    -> suppression threshold   抑制閾值
% N     -> clone number multiplier  群體大小
% Nc    -> no. of clones to be generated  每個候選克隆數(shù)目
% beta  -> decay of the inverse exponential function  延遲轉(zhuǎn)換指數(shù)（衰變的逆指數(shù)函數(shù)）
% gen   -> maximum number of generations   迭代的最大代數(shù)
%

% Function to be Optimized 
clc;
ts = 0.1; N = 20;	Nc = 10; beta = 100;gen = 2000;
f='exp(-2.0*log(2).*(x-0.1).^2/0.64).*(sin(5.*pi*x).^6)';

%f='(sin(5*pi*(x.^0.75-0.05))).^6';
% Parameters for Ploting the Affinity Landscape  親合力參數(shù)
xmin = 0; xmax = 1; x=xmin:0.001:xmax;
vxp = x;
vyp = eval(f);  % 計算親和度

% Initial Random Population Within the Intervals (xmin/xmax; ymin/ymax)
% 間隔之內(nèi)初始化隨機數(shù)
Ab = xmin + rand(N,1).*(xmax - xmin);
x = Ab;
fit = eval(f);  % 計算親和度
figure(1); 
imprime(f,1,vxp,vyp,x,fit,1,1);
title('Initial Population');

% CHECK POSSIBILITY OF USING THE MST TO DETERMINE THE NUMBER OF OPTIMUM ...
   
it = 0; 
Nold = N + 1; 
Nsup = N;
FLAG = 0;   % 找到最優(yōu)解時標志為非零,跳出循環(huán)
FLAGERROR = 0;
avfitold = mean(fit);   % mean()求平均值
avfit = avfitold-1;     
vout = []; 
vavfit = []; 
vN = [];    % 歷代的個體數(shù)目


% Main Loop
while it < gen & FLAG == 0
   % Reproduction (Cloning), Affinity Maturation成熟親和力, and Selection Within Each Clone
     [Ab] = clone_mut_select(Ab,Nc,beta,norma(fit),xmin,xmax,f);     
   % Immune Network Interactions After a Number of Iterations 疊代
   if rem(it,5) == 0
      if abs(1-avfitold/avfit) < .001
         [Ab] = suppress(Ab,ts);    % 抑制操作
         FLAGERROR = 1;
         Nsupold = Nsup;    % 就個體的數(shù)目
         Nsup = size(Ab,1); % 更新個體的數(shù)目
         vN = [vN,Nsup];
         % Convergence Criterion
         if (Nsupold-Nsup) == 0, % & rem(it,20) == 0,
            FLAG = 1; 
            FLAGERROR = 0;
         end;
      end;
   end;
   
   % Insert randomly generated individuals
   if FLAGERROR == 1,
      d = round(.4*N);
      Ab1 = xmin + rand(d,1).*(xmax - xmin);
      Ab = [Ab;Ab1];
      FLAGERROR = 0;
   end;
      
   % Evaluating Fitness
  
   x = Ab;
   fit = eval(f); 
   avfitold = avfit;    % 前一次平均值
   [out,I] = max(fit); 
   avfit = mean(fit);   
   
   % Ploting Results
   imprime(f,1,vxp,vyp,x,fit,it,10);
   N = size(Ab,1);  % 更新N

   it = it + 1; 
   vout = [vout,out]; 
   vavfit = [vavfit,avfit]; % vN = [vN,N];
   disp(sprintf('It: %d	Max: %f	Av: %f	Net size: %d',it,out,avfit,N));
end;



imprime2(f,1,vxp,vyp,x,fit,1,1);
figure(2); plot(vout); hold on; plot(vavfit,'-.'); title('Fitness'); hold off;
figure(3); plot(vN); title('N');

disp(sprintf('找到以下%d個極值點:',N));
ymax=-1000;imax=-1;xmax=-1000;
for i=1:N,
   x = Ab(i); 
   y = eval(f);
   if y>ymax,
       imax=i;
       xmax=x;
       ymax=y;
   end;
disp(sprintf('f(%f)=%f',x,y));
end;

disp(sprintf('其中,最大值f(%f)=%f',xmax,ymax));


%--------------------------------------------------------------------------
% ------------------- %
% SECONDARY FUNCTIONS %
% ------------------- %
%--------------------------------------------------------------------------
function [C] = clone_mut_select(Ab,Nc,beta,fitin,xmin,xmax,f);
%！！！！！！！！！！！！
% 構(gòu)造新的群體
%！！！！！！！！！！！！
% C   -> matrix of clones 克隆矩陣
% g   -> vector with Gaussian mutation) 突變帶菌者
% Ab  -> matrix of antibodies 抗體矩陣
% N   -> cardinality of Ab 秩
% Nc  -> number of clones for each candidate 每個候選克隆數(shù)目
[N,L] = size(Ab);
C = [];
for i = 1:N,
   vones = ones(Nc,1);
   Cc = vones * Ab(i,:);
   g = (randn(Nc,L)./beta) .* exp(-fitin(i)); % 變異程度隨親和度增加而減少
   g(1,:) = zeros(1,L);	% Keep one previous individual for each clone unmutated 保持一個先前完整的個體
   c = Cc + g;  % 制造新的Nc-1個新個體
   % Keeps all elements of the population within the allowed bounds
   % 允許范圍內(nèi)保持克隆的所有元素
   Ixmin = find(c(:,1) < xmin); 
   Ixmax = find(c(:,1) > xmax);
  % Iymin = find(c(:,2) < ymin); Iymax = find(c(:,2) > ymax);
   if ~isempty(Ixmin),
      c(Ixmin,1) = Cc(length(Ixmin),1);
   end;
   if ~isempty(Ixmax),
      c(Ixmax,1) = Cc(length(Ixmax),1);
   end;
   x = c;
   fit = eval(f);
   [out,I] = max(fit);
   C = [C;c(I,:)];  % C contains only the best individuals of each clone 克隆最好的個體
end;

%--------------------------------------------------------------------------
% Function suppress self-recognizing and non-stimulated Ab from Memory (M)
function [M] = suppress(M,ts);
%！！！！！！！！！！！！
% 抑制：計算矩陣M自身的歐拉距離，刪除距離在區(qū)間(0,ts)內(nèi)的個體。
%！！！！！！！！！！！！
% M   -> memory matrix
% D1  -> idiotypic affinity matrix 突變親和力矩陣
D1 = dist(M,M'); % 計算歐拉距離
aux = triu(D1,1); % 取矩陣D1的對角線
[Is,Js] = find(aux>0 & aux<ts);
if ~isempty(Is),
   Is = ver_eq(Is);
   M = extract(M,Is);
   % D1 = extract(D1,Is);
end;

%--------------------------------------------------------------------------
% Search for repeated indexes
function [Is] = ver_eq(I);
%！！！！！！！！！！！！
% 去除數(shù)組I中相同的元素
%！！！！！！！！！！！！
l = length(I); 
Is = [];
if l > 1,
   for i = 1:l-1,
      aux = I(i);
      auxI = I(i+1:end);
      el = find(auxI == aux);
      if isempty(el),
         Is = [Is,aux];
      end;
   end;
   Is = [Is,I(end)];
else,
   Is = I;
end;

%--------------------------------------------------------------------------
% Function Extracts lines from M indexed by I
function [M] = extract(M,I);    % 篩選功能
%！！！！！！！！！！！！
% 從矩陣M中刪除第I行
%！！！！！！！！！！！！
Maux = zeros(size(M));
Maux(I,:) = M(I,:);
M = M - Maux;
[I] = find(M(:,1)~=0);
M = M(I,:);


%--------------------------------------------------------------------------
% Function normalizes matrix over [0,1] 規(guī)格化矩陣
function [Dn] = norma(D);
%！！！！！！！！！！！！
% 把矩陣D里的數(shù)按相應行的最大最小值轉(zhuǎn)化到[0,1]區(qū)間去
%！！！！！！！！！！！！
% Dn  -> normalized vector over [0,1]  規(guī)格化帶菌者
[np,ni] = size(D);
if ni == 1,
   Dn = (D - min(D))./(max(D)-min(D));  
else,
   vmaxD = max(D); 
   vminD = min(D);
   for i=1:ni,
      Dn(:,i) = (D(:,i) - vminD(i))./(vmaxD(i)-vminD(i));
   end;
end;
% End Function NORMA

% Print Affinity Landscape and Population of Individuals 畫出親和力以及個體數(shù)目

%--------------------------------------------------------------------------
function [] = imprime(f,PRINT,vx,vy,x,fx,it,mit);
%！！！！！！！！！！！！
% 繪制x-fx在f上的分布
%！！！！！！！！！！！！
% x,fx				-> current values
% vxplot, vplot	-> original (base) function
if PRINT == 1,
   if rem(it,mit) == 0, % rem()為求it除mit的余數(shù)
      hold on; 
      fplot(f,[0 1 0 1],'k-')
      xlabel('x');  ylabel('f(x)');
      plot(x,fx,'k*'); drawnow; hold off; % pause
   end;
end;

%--------------------------------------------------------------------------
function [] = imprime2(f,PRINT,vx,vy,x,fx,it,mit);
%！！！！！！！！！！！！
% 繪制x-fx在f上的分布
%！！！！！！！！！！！！
% x,fx				-> current values
% vxplot, vplot	-> original (base) function
if PRINT == 1,
   if rem(it,mit) == 0,
      hold on; 
      fplot(f,[0 1 0 1],'k-')
      xlabel('x');  ylabel('f(x)');
      plot(x,fx,'ro'); drawnow; hold off; % pause
   end;
end;