function [bestnest,fmin]=cuckoo_search(n)

global centroids

Lb = [1 30];
Ub = [100 40];

if nargin<1,
% Number of nests (or different solutions)
n = 25;
end

% Discovery rate of alien eggs/solutions
pa = 0.25;
%D;T1;T2;W1;W2

%% Change this if you want to get better results
% Tolerance
%Tol = 3;

% Número máximo de iterações
Nmi = 100;
%% Simple bounds of the search domain

% Random initial solutions
for i=1:n
    nest(i,:)=Lb+(Ub-Lb).*rand(size(Lb));
end

% Get the current best
fitness=10^10*ones(n,1);
[fmin,bestnest,nest,fitness]=get_best_nest(nest,nest,fitness);
N_iter=0;
%% Starting iterations
c = 1;
while (c < Nmi+1)
    % Generate new solutions (but keep the current best)
     new_nest=get_cuckoos(nest,bestnest,Lb,Ub);   
     [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
     
    % Update the counter
      N_iter=N_iter+n; 
      
    % Discovery and randomization
      new_nest=empty_nests(nest,Lb,Ub,pa);
    
    % Evaluate this set of solutions
      [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
      
    % Update the counter again
      N_iter=N_iter+n;
      
    % Find the best objective so far  
    if fnew<fmin,
        fmin=fnew;
        bestnest=best
        save('C:\Users\flavi\Documents\DOUTORADO\5G Cluster\gabriel\saved data\BestCentroidsCS3500.txt', 'centroids', '-ASCII');
    end
    
    vetorbest(c) = fmin;
    vetormediafmin(c) = mean(fitness);
    %save('C:\Users\flavi\Documents\MATLAB\Comp Urbana\VetorBestCucoLora.txt', 'vetorbest', '-ASCII');
    %save('C:\Users\flavi\Documents\MATLAB\Comp Urbana\VetorMediaCucoLora.txt', 'vetormediafmin', '-ASCII');

    c = c + 1;
end %% End of iterations

     plot(1:c-1,vetorbest,'-b','LineWidth',2); 
     hold on
     plot(1:c-1,vetormediafmin,'k-','LineWidth',2);
     hold on
     fminmedio=0;
     legend('Best Solution','Mean Solution');
     xlabel('Iteration');
     ylabel('Fitness');

%plot(1:c-1,vetorbest,'-b',1:c-1,vetormediafmin,'k-');
%legend({'Best Fitness Value','Average Fitness Value'}, 'FontSize', 12) 

%% Post-optimization processing
%% Display all the nests
% disp(strcat('Total number of iterations=',num2str(N_iter)));
fmin

%%%%%%%% PLOT %%%%%%%%%%

bestnest
%save('C:\Users\flavi\Documents\MATLAB\Comp Urbana\BestLora.txt', 'bestnest', '-ASCII');


%% --------------- All subfunctions are list below ------------------
%% Get cuckoos by ramdom walk
function nest=get_cuckoos(nest,best,Lb,Ub)
% Levy flights
n=size(nest,1);

% Levy exponent and coefficient
% For details, see equation (2.21), Page 16 (chapter 2) of the book
% X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).

beta=1/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);

for j=1:n,
    s=nest(j,:);
    % This is a simple way of implementing Levy flights
    % For standard random walks, use step=1;
    %% Levy flights by Mantegna's algorithm
    u=randn(size(s))*sigma;
    v=randn(size(s));
    step=u./abs(v).^(1/beta);
  
    % In the next equation, the difference factor (s-best) means that 
    % when the solution is the best solution, it remains unchanged.     
    stepsize=0.01*step.*(s-best);
    % Here the factor 0.01 comes from the fact that L/100 should the typical
    % step size of walks/flights where L is the typical lenghtscale; 
    % otherwise, Levy flights may become too aggresive/efficient, 
    % which makes new solutions (even) jump out side of the design domain 
    % (and thus wasting evaluations).
    % Now the actual random walks or flights
    s=s+stepsize.*randn(size(s));
   % Apply simple bounds/limits
   nest(j,:)=simplebounds(s,Lb,Ub);
end

end

%% Find the current best nest
function [fmin,best,nest,fitness]=get_best_nest(nest,newnest,fitness)
% Evaluating all new solutions
for j=1:size(nest,1),
    fnew=fobj(newnest(j,:));
    if fnew<=fitness(j),
       fitness(j)=fnew;
       nest(j,:)=newnest(j,:);
    end
end
% Find the current best
[fmin,K]=min(fitness);
best=nest(K,:);

end

%% Replace some nests by constructing new solutions/nests
function new_nest=empty_nests(nest,Lb,Ub,pa)
% A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
% Discovered or not -- a status vector
K=rand(size(nest))>pa;

% In the real world, if a cuckoo's egg is very similar to a host's eggs, then 
% this cuckoo's egg is less likely to be discovered, thus the fitness should 
% be related to the difference in solutions.  Therefore, it is a good idea 
% to do a random walk in a biased way with some random step sizes.  
%% New solution by biased/selective random walks
stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));
new_nest=int16(nest+stepsize.*K);
for j=1:size(new_nest,1)
    s=new_nest(j,:);
  new_nest(j,:)=simplebounds(s,Lb,Ub);  
end

end

% Application of simple constraints
function s=simplebounds(s,Lb,Ub)
  % Apply the lower bound
  ns_tmp=s;
  I=ns_tmp<Lb;
  ns_tmp(I)=Lb(I);
  
  % Apply the upper bounds 
  J=ns_tmp>Ub;
  ns_tmp(J)=Ub(J);
  % Update this new move 
  s=ns_tmp;
  
end

%% You can replace the following by your own functions
% A d-dimensional objective function
function z=fobj(u)

global populacao_data2;

w=[0.7 0.3];

inicio=int16(u(1));
    
dados22(inicio);

for i=1:length(centroids)
    numerocentroide=i;  %para o python deveria ser i-1
    num=1;
    contador=0;
    for j=1:length(populacao_data2),
        if numerocentroide==populacao_data2(j,3),
            if ((centroids(i,1)~=populacao_data2(j,2)) && (centroids(i,2)~= populacao_data2(j,1))),
            hav=haversine([centroids(i,1) populacao_data2(j,2)],[centroids(i,2) populacao_data2(j,1)]);

            distanciahaversine=hav;%metros
            
            x=[13.958 5.8 0.759 1.862]; 
            ht=40;%40;%metros
            hr=1.6;%%1.6;%metros
            f=3.5e9;%Hertz
            pl=ecc33(x,ht,hr,f,distanciahaversine);
            
            pt=u(2);%dBm
            gt=3;%15.6;%dBi
            gr=0; %dBi
            pr(num)=pt-pl+gt+gr;%dBm
          
            
            if pr(num)>=-90 
                contador=contador+1;
            end; 
            
            num=num+1;
            
            end
        end
    end
    %percentualdecoberturaporcentroide(i)=(contador/length(pr))*100;
    percentualdecoberturaporcentroide(i)=(contador/num)*100;
%     mediapotencia(i)=mean(pr);
    clear pr
end

   mediadopercentualdecoberturadetodososclusters=mean(percentualdecoberturaporcentroide);
   %clear mediapotencia;
   clear percentualdecoberturaporcentroide;

    centroids;

    mediadopercentualdecoberturadetodososclusters;

    obj1 = 100 - mediadopercentualdecoberturadetodososclusters;
    obj2 = abs(pt-30);

    z=(obj1*w(1,1)+obj2*w(1,2))

end

end