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

% function [M] = optd2(ts,N,Nc,beta,gen)
% M     -> matrix of memory cells 記憶細胞矩陣
% ts    -> suppression threshold   抑制極限
% N     -> clone number multiplier  克隆的數目系數
% Nc    -> no. of clones to be generated  克隆所產生的個體數 
% beta  -> decay of the inverse exponential function  延遲轉換指數
% gen   -> maximum number of generations  生成的最大數目
%

% Function to be Optimized 
clc;
ts = 0.1; N = 20;	Nc = 10; beta = 100;gen = 2000;
%f = '1 * x .* sin(4 * pi .* x) - 1 * y.* sin(4 * pi .* y + pi) + 1'; 
%xmin = -1; xmax = 2; ymin = -1; ymax = 2;
%[x,y] = meshgrid(xmin:0.03:xmax,ymin:0.03:ymax);

f='100. *(x.^2 - y).^2+(1-x).^2';
xmin = -2.048; xmax = 2.048; ymin = -2.048; ymax = 2.048;
[x,y] = meshgrid(xmin:0.01:xmax,ymin:0.01:ymax); 

vxp = x; vyp = y;
vzp = eval(f);

% Initial Random Population Within the Intervals (xmin/xmax; ymin/ymax)
Ab1 = xmin + rand(N,1).*(xmax - xmin);
Ab2 = ymin + rand(N,1).*(ymax - ymin);
Ab = [Ab1,Ab2];
x = Ab(:,1); y = Ab(:,2);
fit = eval(f);
figure(1); imprime(1,vxp,vyp,vzp,x,y,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; FLAGERROR = 0;
avfitold = mean(fit); avfit = avfitold-1;
vout = []; vavfit = []; vN = [];
% 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,ymin,ymax,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; Nsup = size(Ab,1); 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);
      Ab2 = ymin + rand(d,1).*(ymax - ymin);
      Ab = [Ab;Ab1,Ab2];
      FLAGERROR = 0;
   end;
      
   % Evaluating Fitness
   x = Ab(:,1); y = Ab(:,2);
   z = complex(x,y);
   fit = eval(f); avfitold = avfit;
   [out,I] = max(fit); avfit = mean(fit);
   
   % Ploting Results
   imprime(1,vxp,vyp,vzp,x,y,fit,it,10);
   N = size(Ab,1);

   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;
imprime(1,vxp,vyp,vzp,x,y,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));
zmax=-1000;imax=-1;xmax=-1000;ymax=-1000;
for i=1:N,
   x = Ab(i,1); y = Ab(i,2);
   z = eval(f);
   if z>zmax,
       imax=i;xmax=x;ymax=y;zmax=z;
   end;
disp(sprintf('f(%f,%f)=%f',x,y,z));
end;

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


% ------------------- %
% SECONDARY FUNCTIONS %
% ------------------- %

function [C] = clone_mut_select(Ab,Nc,beta,fitin,xmin,xmax,ymin,ymax,f);
% C   -> matrix of clones
% g   -> vector with Gaussian mutation
% Ab  -> matrix of antibodies
% N   -> cardinality of Ab
% Nc  -> number of clones for each candidate
[N,L] = size(Ab);
C = [];
for i=1:N,
   vones = ones(Nc,1);
   Cc = vones * Ab(i,:);
   rrr=randn(Nc,L);
    g = (rrr./beta) .* exp(-fitin(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;
   % Keeps all elements of the population within the allowed bounds
   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;
   if ~isempty(Iymin),
      c(Iymin,2) = Cc(length(Iymin),2);
   end;
   if ~isempty(Iymax),
      c(Iymax,2) = Cc(length(Iymax),2);
   end;
   x = c(:,1); y = c(:,2);
   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   -> memory matrix
% D1  -> idiotypic affinity matrix
D1 = dist(M,M');
aux = triu(D1,1);
[Is,Js] = find(aux>0 & aux<ts);
if ~isempty(Is),
   Is = ver_eq(Is);
   M = extract(M,Is);
   % D1 = extract(D1,Is);
end;
% D1 = dist(M,M');

% Search for repeated indexes
function [Is] = ver_eq(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);
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]
function [Dn] = norma(D);
% Dn  -> normalized vector over [0,1]
[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
function [] = imprime(PRINT,vx,vy,vz,x,y,fx,it,mit);
% x,fx				-> current values
% vxplot, vplot	-> original (base) function
if PRINT == 1,
   if rem(it,mit) == 0,
     mesh(vx,vy,vz); 
      hold on; %axis([-1 2 -1 2 -2 4]);
      xlabel('x'); ylabel('y'); zlabel('f(x,y)');
      plot3(x,y,fx,'k*'); drawnow; hold off; % pause
   end;
end;