% Firefly Algorithm by X-S Yang (Cambridge University)
% Multidimensional implementation prepared by Szymon Lukasik and Slawomir Zak
% Bibliographical reference: S. Lukasik, S. Zak, "Firefly Algorithm for Continuous Constrained Optimization", Lecture Notes in Artificial Intelligence, vol. 5796, pp. 97-106, 2009.

function [stats]=firefly(seed,visual,eval_function,global_opt,global_error,n,steps,alpha,gamma,limits, beta0)
% 'eval_function' string of eval function
% n - number of fireflies
% steps - number of steps
% alpha - randomness [0,1]
% gamma - absorption coefficient
% [limits] - matrix determining variables limits

%default values
def_rand = randseed;
def_visual = true;
def_eval_function = 'four_peaks';
def_global_opt = -2.0;
def_n=12;
def_steps = 50;
def_alpha = 0.05;
def_gamma = 1.0;
def_limits = [-5 5 -5 5];
def_error = 1e-6;
def_beta0 = 1;

switch nargin
    case 0
        seed=def_rand; visual=def_visual; eval_function=def_eval_function; global_opt=def_global_opt; global_error=def_error; n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 1
        visual=def_visual; eval_function=def_eval_function; global_opt=def_global_opt; global_error=def_error; n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0; 
    case 2
        eval_function=def_eval_function; global_opt=def_global_opt; global_error=def_error; n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits;     
    case 3
        global_opt=def_global_opt; global_error=def_error; n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 4
        global_error=def_error; n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 5 
        n=def_n; steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 6 
        steps=def_steps; alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 7 
        alpha=def_alpha; gamma=def_gamma; limits = def_limits; beta0 = def_beta0;
    case 8 
        gamma=def_gamma; beta0 = def_beta0;
    case 9 
        limits = def_limits; beta0 = def_beta0;   
    case 10
        beta0 = def_beta0;
        
        
end
best_qual=inf;

rand('state', seed);
% plot if two-dimensional problem and data visualization enabled
if ((length(limits)/2==2) && (visual == true))     
    grid_size=100;
    dx=(limits(2)-limits(1))/grid_size;
    dy=(limits(4)-limits(3))/grid_size;
    [x_grid,y_grid]=meshgrid(limits(1):dx:limits(2),limits(3):dy:limits(4));    
    z=zeros(length(limits(1):dx:limits(2)),length(limits(3):dx:limits(4)));
    for i=1:length(limits(1):dx:limits(2))
        for j=1:length(limits(3):dx:limits(4))
            z(i,j)=feval(eval_function,[x_grid(i,j),y_grid(i,j)]);
        end
    end
    figure(1); surfc(x_grid,y_grid,z); 
    figure(2);
    contour(x_grid,y_grid,z);  
end
% end of contour plotting

% generate initial solution
[x_new,x_qual_new]=Generate_Initial_Sol(n,limits);
iter = 1; optym=false;
startTime = cputime;
for i=1:steps
  for j=1:n
      x_qual_new(j)=feval(eval_function,x_new(j,:));
  end
  % sort fireflies
  [x_qual_new,light_index]=sort(x_qual_new);
  x_new=x_new(light_index,:);
  if best_qual>x_qual_new(1)
      best=x_new(1,:);
      best_qual=x_qual_new(1);
  end
  x_old=x_new;
  x_qual_old=x_qual_new;
  % calculate new configuration
  [x_new]=Move_Fireflies(x_new,x_qual_new, x_old, x_qual_old, alpha, gamma, limits, beta0);      
  if (abs(x_qual_new(1)-global_opt) <= global_error), optym=true; break; end;
  
  % plot if two-dimensional problem and data visualization enabled
  if ((length(limits)/2==2) && (visual == true))    
      contour(x_grid,y_grid,z);
      hold on
      plot(x_new(:,1),x_new(:,2),'.');
      axis (limits);
      M(i) = getframe;
      hold off
  end;
  iter = iter + 1;
end
endTime = cputime - startTime;

stats = [best best_qual optym iter-1 endTime];


function [x_new,x_qual_new]=Generate_Initial_Sol(n,limits)
  x_new=zeros(n,length(limits)/2);
   for i=1:length(limits)/2
    x_limit=limits(2*i)-limits(2*i-1);
    x_new(:,i)=rand(n,1)*x_limit+limits(2*i-1);
  end
  x_qual_new=zeros(n,1);
  
function [x_new]=Move_Fireflies(x_new,x_qual_new, x_old, x_qual_old, alpha, gamma, limits, beta0)
pop=size(x_old,1);
dimension=size(x_old,2);

for i=1:pop
  for j=1:pop
    r=norm(x_new(i,:)-x_old(j,:));
    if x_qual_new(i)>x_qual_old(j)        
        beta=beta0*exp(-gamma*r^2);
        % generate random vector (randomization factor)
        random_vector=rand(1,dimension);
        for k=1:dimension
            if sign(rand(1,1) - 0.5)<0
                random_vector(k)=random_vector(k)*(limits(2*k)-x_new(i,k));
            else
                random_vector(k)=-random_vector(k)*(x_new(i,k)-limits(2*k-1));
            end
        end
        x_new(i,:)=x_new(i,:)*(1-beta)+x_old(j,:)*beta+alpha.*random_vector;
    end
  end
end

% and now move randomly the best firefly
if ispref('firefly','move_best')==1 && getpref('firefly','move_best')==true 
    best_index = 1;
    % generate random vector (randomization factor)
    random_vector=rand(1,dimension);
    for k=1:dimension
     %if sign(random_vector(k)-0.5)<0
     if sign(rand(1,1)-0.5)<0
         random_vector(k)=random_vector(k)*(limits(2*k)-x_new(best_index,k));
     else
         random_vector(k)=-random_vector(k)*(x_new(best_index,k)-limits(2*k-1));
     end
    end
    x_new(best_index,:)=x_new(best_index,:)+alpha.*random_vector;
end


function y=four_peaks(x)
y=-(exp(-(x(1)-4)^2-(x(2)-4)^2)+exp(-(x(1)+4)^2-(x(2)-4)^2)+2*exp(-x(1)^2-(x(2)+4)^2)+2*exp(-x(1)^2-x(2)^2));


function y=himmelblau(x)  % [-6 6 -6 6] f(x)=0
y=(x(1)^2+x(2)-11)^2 + (x(1)+x(2)^2-7)^2; 


function y=beale(x) % [?4.5 4.5 -4.5 4.5]  The global minimum: x = (3, 0.5), f(x) = 0.
y=(1.5-x(1)*(1-x(2)))^2+(2.25-x(1)*(1-x(2)^2))^2+(2.625-x(1)*(1-x(2)^3))^2;

% ZMIANA z 2 na 10
function y = ackley(x)
n = 10; 
a = 20; b = 0.2; c = 2*pi;
s1 = 0; s2 = 0;
for i=1:n;
   s1 = s1+x(i)^2;
   s2 = s2+cos(c*x(i));
end
y = -a*exp(-b*sqrt(1/n*s1))-exp(1/n*s2)+a+exp(1);

% ZMIANA z 2 na 10
function y = schwel(x) % [-500 500 -500 500] f(1,1) = 0;
n = 2; 
s = sum(-x.*sin(sqrt(abs(x))));
y = 418.9829*n+s;


function y = griewank(x)
n = 5;
fr = 4000;
s = 0;
p = 1;
for j = 1:n; s = s+x(j)^2; end
for j = 1:n; p = p*cos(x(j)/sqrt(j)); end
y = s/fr-p+1;


function y = rastrigin(x) % The global minima: x* =  (0, …, 0), f(x*) = 0 domain: [-5.12 5.12 -5.12 5.12]
n = 20; 
s = 0;  
for j = 1:n
    s = s+(x(j)^2-10*cos(2*pi*x(j))); 
end
y = 10*n+s;

% ZMIANA z 2 na 4
function y = powersum(x)
n = 4;  
b = [8,18,44,114];
s_out = 0;
for k = 1:n;
    s_in = 0;
    for j = 1:n
        s_in = s_in+x(j)^k;
    end
    s_out = s_out+(s_in-b(k))^2;
end
y = s_out;

% ZMIANA z 2 na 5
function y = michalewicz(x) % [ 0 pi 0 pi ] 
n = 2;  
m = 10;
s = 0;
for i = 1:n;
    s = s+sin(x(i))*(sin(i*x(i)^2/pi))^(2*m);
end
y = -s;

% ZMIANA z 2 na 3
function y = sphere(x) % [-5.12 5.12 -5.12 5.12]
n = 2; 
s = 0;
for j = 1:n
    s = s+x(j)^2; 
end
y = s;

% ZMIANA z 2 na 5
function y = levy(x) % [ -10 10 -10 10]
n = 2;
for i = 1:n; z(i) = 1+(x(i)-1)/4; end
s = sin(pi*z(1))^2;
for i = 1:n-1
    s = s+(z(i)-1)^2*(1+10*(sin(pi*z(i)+1))^2);
end 
y = s+(z(n)-1)^2*(1+(sin(2*pi*z(n)))^2);


function y = branin(x)  %[-5 10 -5 10] f(x*) = 0.397887.
y = (x(2)-(5.1/(4*pi^2))*x(1)^2+5*x(1)/pi-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;


function y = shekel(x) %[0 10 0 10 0 10 0 10]
m = 10; % number of minima
a = ones(10,4);
a(1,:) = 4.0*a(1,:);
a(2,:) = 1.0*a(2,:);
a(3,:) = 8.0*a(3,:);
a(4,:) = 6.0*a(4,:);
for j = 1:2;
   a(5,2*j-1) = 3.0; a(5,2*j) = 7.0; 
   a(6,2*j-1) = 2.0; a(6,2*j) = 9.0; 
   a(7,j)     = 5.0; a(7,j+2) = 3.0;
   a(8,2*j-1) = 8.0; a(8,2*j) = 1.0;
   a(9,2*j-1) = 6.0; a(9,2*j) = 2.0;
   a(10,2*j-1)= 7.0; a(10,2*j)= 3.6;
end
c(1) = 0.1; c(2) = 0.2; c(3) = 0.2; c(4) = 0.4; c(5) = 0.4;
c(6) = 0.6; c(7) = 0.3; c(8) = 0.7; c(9) = 0.5; c(10)= 0.5;
s = 0;
for j = 1:m;
   p = 0;
   for i = 1:4
      p = p+(x(i)-a(j,i))^2;
   end
   s = s+1/(p+c(j));
end
y = -s;

function y = hartman6(x) %[0 1 0 1 0 1 0 1 0 1 0 1 ]
a(1,1)=10.0;	a(1,2)=3.0;		a(1,3)=17.0;	a(1,4)=3.5;		a(1,5)=1.7;		a(1,6)=8.0;
a(2,1)=0.05;	a(2,2)=10.0;	a(2,3)=17.0;	a(2,4)=0.1;		a(2,5)=8.0;		a(2,6)=14.0;
a(3,1)=3.0;		a(3,2)=3.5;		a(3,3)=1.7;		a(3,4)=10.0;	a(3,5)=17.0;	a(3,6)=8.0;
a(4,1)=17.0;	a(4,2)=8.0;		a(4,3)=0.05;	a(4,4)=10.0;	a(4,5)=0.1;		a(4,6)=14.0;
c(1)=1.0;c(2)=1.2;c(3)=3.0;c(4)=3.2;
p(1,1)=0.1312;	p(1,2)=0.1696;	p(1,3)=0.5569;	p(1,4)=0.0124;	p(1,5)=0.8283;	p(1,6)=0.5886;
p(2,1)=0.2329;	p(2,2)=0.4135;	p(2,3)=0.8307;	p(2,4)=0.3736;	p(2,5)=0.1004;	p(2,6)=0.9991;
p(3,1)=0.2348;	p(3,2)=0.1451;	p(3,3)=0.3522;	p(3,4)=0.2883;	p(3,5)=0.3047;	p(3,6)=0.6650;
p(4,1)=0.4047;	p(4,2)=0.8828;	p(4,3)=0.8732;	p(4,4)=0.5743;	p(4,5)=0.1091;	p(4,6)=0.0381;
s = 0;
for i=1:4;
   sm=0;
   for j=1:6;
      sm=sm+a(i,j)*(x(j)-p(i,j))^2;
   end
   s=s+c(i)*exp(-sm);
end
y = -s;

%mam
function y = hartman3(x) %[0 1 0 1 0 1]
a(:,2)=10.0*ones(4,1);
for j=1:2;
   a(2*j-1,1)=3.0; a(2*j,1)=0.1; 
   a(2*j-1,3)=30.0; a(2*j,3)=35.0; 
end
c(1)=1.0;c(2)=1.2;c(3)=3.0;c(4)=3.2;
p(1,1)=0.36890;p(1,2)=0.11700;p(1,3)=0.26730;
p(2,1)=0.46990;p(2,2)=0.43870;p(2,3)=0.74700;
p(3,1)=0.10910;p(3,2)=0.87320;p(3,3)=0.55470;
p(4,1)=0.03815;p(4,2)=0.57430;p(4,3)=0.88280;
s = 0;
for i=1:4;
   sm=0;
   for j=1:3;
      sm=sm+a(i,j)*(x(j)-p(i,j))^2;
   end
   s=s+c(i)*exp(-sm);
end
y = -s;

%mam
function y = permutation(x) %[-4 4 -4 4]
n = 2;
b = 10;
s_out = 0;
for k = 1:n;
    s_in = 0;
    for j = 1:n
        s_in = s_in+(j+b)*(x(j)^k-(1/j)^k);
    end
    s_out = s_out+s_in^2;
end
y = s_out;

% ZMIANA z 2 na 4
%mam
function y = rosenbrock(x) % [ -5 10 -5 10]
n = 4;
sum = 0;
for j = 1:n-1;
    sum = sum+100*(x(j)^2-x(j+1))^2+(x(j)-1)^2;
end
y = sum;

%mam
function y = easom(x) %[-100 100 -100 100] f(x*) = -1
y = -cos(x(1))*cos(x(2))*exp(-(x(1)-pi)^2-(x(2)-pi)^2);