function [best,fmin,N_iter]=fpa_demo(para)
clear all
close all
clc
global w;

global centroids;

% Default parameters
if nargin<1,
   para=[25 0.50];
end

n=para(1);           % Population size, typically 10 to 25
p=para(2);           % probabibility switch

% Iteration parameters
N_iter= 100 ;            % Total number of iterations %2000

% Dimension of the search variables

tic;

NPareto=2;  % Number of points on the Pareto front 

for kpareto=1:NPareto,
    % Generate a weighting coefficient:w so that w1=w, w2=1-w, w1+w2=1.  
    % Observations suggest that systematically monotonic weights are
    % better than random weights.
    w=[0.7 0.3]; 


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

d=length(Lb);
for i=1:n,
  Sol(i,:)=Lb+(Ub-Lb).*rand(1,d);
  Fitness(i)=Fun(Sol(i,:)); 
end

% Find the current best
[fmin,I]=min(Fitness);
best=Sol(I,:);
S=Sol; 
fminmedio=0;

% Start the iterations -- Flower Algorithm 
for t=1:N_iter,
        % Loop over all bats/solutions
        for i=1:n,
          % Pollens are carried by insects and thus can move in
          % large scale, large distance.
          % This L should replace by Levy flights  
          % Formula: x_i^{t+1}=x_i^t+ L (x_i^t-gbest)
          if rand>p,
          %% L=rand;
          L=Levy(d);
          
%           if t==1
%               Lvetor(i,:)=L/0.01;
%               figure(6)
%               hold on
%               plot(Lvetor(:,1),'*-k','LineWidth',2)
%               plot(Lvetor(:,2),'*-b','LineWidth',2)
%           end;
          dS=L.*(Sol(i,:)-best);
          S(i,:)=Sol(i,:)+dS;
          
          % Check if the simple limits/bounds are OK
          S(i,:)=simplebounds(S(i,:),Lb,Ub);
          
          % If not, then local pollenation of neighbor flowers 
          else
              epsilon=rand;
              % Find random flowers in the neighbourhood
              JK=randperm(n);
              % As they are random, the first two entries also random
              % If the flower are the same or similar species, then
              % they can be pollenated, otherwise, no action.
              % Formula: x_i^{t+1}+epsilon*(x_j^t-x_k^t)
              S(i,:)=S(i,:)+epsilon*(Sol(JK(1),:)-Sol(JK(2),:));
              % Check if the simple limits/bounds are OK
              S(i,:)=simplebounds(S(i,:),Lb,Ub);
          end
          
          % Evaluate new solutions
           Fnew=Fun(S(i,:));
          % If fitness improves (better solutions found), update then
            if (Fnew<=Fitness(i)),
                Sol(i,:)=S(i,:);
                Fitness(i)=Fnew;
           end
           
          % Update the current global best
          if Fnew<=fmin,
                best=S(i,:);
                fmin=Fnew;
          save('C:\Users\flavi\Documents\DOUTORADO\5G Cluster\gabriel\saved data\BestCentroidsFPA3500.txt', 'centroids', '-ASCII');
          end
          
          fminmedio=Fnew+fminmedio;
        end
        % Display results every 100 iterations
        if round(t/100)==t/100,
        best;
        fmin;
        end
        
    v_fminmedio(t)=fminmedio/n;    
    v_best1(t)=fmin;
    %w_best(t)=best(:,1);
    %t_best(t)=best(:,2);
    vetorgerac(t)=t;
    
    %save('C:\Users\flavi\Documents\MATLAB\Comp Urbana\FPA\VetorBestFPALora.txt', 'v_best1', '-ASCII');
    %save('C:\Users\flavi\Documents\MATLAB\Comp Urbana\FPA\VetorMediaFPALora.txt', 'v_fminmedio', '-ASCII');
    %plot de todos os melhores valores por geração 
    if kpareto==1
    
     
     plot(vetorgerac,v_best1,'-b','LineWidth',2) 
     hold on
     plot(vetorgerac,v_fminmedio,'-k','LineWidth',2) 
     hold on
     fminmedio=0;
     hi=legend('Best Solution','Mean Solution');
     xlabel('Iteration');
     ylabel('Fitness');
     
    end;
        
end
          
% kpareto
% 
% [v_bestbest,tt]=min(v_best1);
%  v_mediabest(kpareto)=mean(v_best1);
%     v_best2(kpareto)=v_bestbest;
%     w_bestbest(kpareto)=w_best(tt);
%     t_bestbest(kpareto)=t_best(tt);
    
% Output/display
disp(['Total number of evaluations: ',num2str(N_iter*n)]);
disp(['Best solution= ',num2str(int16(best(1,1))),'   ', num2str(best(1,2)),'   fmin= ',num2str(fmin)]);
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;


% Draw n Levy flight sample
function L=Levy(d)
% Levy exponent and coefficient
% For details, see Chapter 11 of the following book:
% Xin-She Yang, Nature-Inspired Optimization Algorithms, Elsevier, (2014).
beta=1/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
    u=randn(1,d)*sigma;
    v=randn(1,d);
    step=u./abs(v).^(1/beta);
    L=0.01*step; 

% Objective function and here we used Rosenbrock's 3D function
function z=Fun(u)
global w;

global centroids populacao_data2;

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 = abs(100 - mediadopercentualdecoberturadetodososclusters);
obj2 = abs(pt-30);

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