clc;
clear all;
close all;

total_nodes = 100;           
node_ids = 1:total_nodes;   


init_thres = 100;   
init_flag = 0;           
end_flg = 0;              
iterations_start = 0;                
while (init_flag == 0)
    [G, init_flag] = createConnectedNetwork (total_nodes);
    iterations_start = iterations_start + 1;    
    if ((iterations_start >= 100) && (init_flag == 0))     
        end_flg = 1;        
        break; 
    end 
end
if (end_flg == 1) 
    disp ('Cannot estabish a network');
    disp ('Try to change values for initiating the network');
    disp ('Network Design stoppped');
else
    
   
    disp (strcat('A connected network was created with number of iterations_start = ', int2str(iterations_start)));

    
    mpr_selector = datasample(node_ids, 1);
    disp (strcat('MPR selector node: ', int2str(mpr_selector)));

    
    [firstHop_ids, secondHop_ids] = get_First_Second_Neighbors(mpr_selector, G);
    disp ('Neighbourhood in First iterative hop:');
    disp (firstHop_ids');
    disp ('Second-hop neighbors:');
    disp (secondHop_ids);
    
    
    if (isempty(secondHop_ids))
        disp ('No second-hop neighbors exist');
        disp('try re-creating the network');
        disp('Program stopped');
    else
        
        
        selected_MPRs = [];
     
        
        for i=1:length(secondHop_ids)  
            
            [included, count_included] = getIncludedNeighbors(secondHop_ids(i), G, firstHop_ids);
            if ((count_included == 1) && (~ismember(included(1), selected_MPRs)))
                
                selected_MPRs = [selected_MPRs, included];
            end      
        end 
       
      
        available_firstHop = firstHop_ids;
        if (~isempty(selected_MPRs))
            disp ('selected MPRs in step 1:');
            disp (selected_MPRs);
           
            selected_ids= ismember(available_firstHop, selected_MPRs);
            available_firstHop(selected_ids)=[];
        else
            disp ('No MPRs selected in step 1');
        end 
        
      
        if (isempty(available_firstHop))
            disp('All first-hop neighbors were selected in step 1');
            disp('Program stopped');
        else 
            disp ('First-hop neighbors available for selection:');
            disp (available_firstHop');
            
            
            
            
            uncovered_secondHop = secondHop_ids;
            
            if (~isempty(selected_MPRs))
                for i=1:length(selected_MPRs)
                    m = selected_MPRs(i);
                   
                    [temp, count_temp] = getIncludedNeighbors(m, G, uncovered_secondHop);
                   
                    ids= ismember(uncovered_secondHop, temp);
                    uncovered_secondHop(ids)=[];
                end 
            end 
        
            if (isempty(uncovered_secondHop))
                disp('All second-hop neighbors are covered');
                disp('Program stopped');
            else
                disp('uncovered second-hop neighbors');
                disp(uncovered_secondHop);
                
                
                while (~isempty(uncovered_secondHop))
                    
                    disp('selecting the first-hop neighbor with max coverage');
                    
                    
                    [selected_node, covered_nodes] = getNodeMaxCoverage(available_firstHop, uncovered_secondHop, G);
                    
                   
                    selected_MPRs = [selected_MPRs, selected_node];
                    
                    disp ('selected MPRs changed to:');
                    disp (selected_MPRs);
                    
               
                    idx= ismember(available_firstHop, selected_MPRs);
                    available_firstHop(idx)=[];
                    
                   
                    ids= ismember(uncovered_secondHop, covered_nodes);
                    uncovered_secondHop(ids)=[];
                    disp ('Uncovered second-hop neighbors changed to:');
                    disp (uncovered_secondHop);
                
                end 
            end 
       
        end 
        
        
        disp('Overall Selected MPRs:');
        disp(selected_MPRs);
        
        
        saving = length(firstHop_ids) - length(selected_MPRs);
        disp('First-hop nodes used for broadcasting decreased by:');
        disp(saving);
        disp ('end');
    
         
        title_s = 'MPR Selector (Red), First-hop (Green), Second-hop (Black), selected MPRs (Yellow)';
    
     
        figure('units','normalized','outerposition',[0 0 1 1])
       
        h = plot(G);
         
        highlight(h, mpr_selector, 'NodeColor', 'r');
        
        highlight(h, firstHop_ids,'NodeColor', 'g');
        
        highlight(h, secondHop_ids,'NodeColor', 'k');
        
        highlight(h, selected_MPRs,'NodeColor', 'y');
        
        title(title_s);
    end  
end 
%%
source_in.x=0.5*total_nodes;  %location of source_in on x-axis
source_in.y=0.5*total_nodes;  %location of source_in on y-axis
n=total_nodes; %nodes
P=0.1;  %probability of cluster heads
Eo=0.5;%initial energy
%
Energy_check_init=Eo;
%
energy_tx_init=50*0.000000001;  %tx energy
energy_rx_init=50*0.000000001;  %rx energy
Efs=10*0.000000000001;  %free space loss
Emp=0.0013*0.000000000001;   %multipath loss
%Data Aggregation Energy
EDA=5*0.000000001;  %compression energy
a=1;   %fraction of energy enhancment of advance nodes
rmax=5000;  %maximum number of rounds
do=hybrid_dist(Efs,Emp);  %distance do is measured based on manhattan and euclidian
Et=0;  %variable just use below 
m=0;
mony=0;
c=0.02;
b=0.7;
talhar=0;
for i=1:1:n
    S(i).xd=covered_nodes(1,1); 
    XR(i)=S(i).xd;
    S(i).yd=covered_nodes(2,1);  
    YR(i)=S(i).yd;
    S(i).G=0; 
    talhar=abs(covered_nodes(1,1)^2*a*10^-4-covered_nodes(1,1)*a*10^-2-covered_nodes(1,1)*10^-2).*(1/sqrt(1-exp(-i)));% energy with modified sigmoid function
    S(i).E=Eo*(1+talhar);
    S(i).A=talhar;
    E(i)= S(i).E;
    if (E(i)>Energy_check_init)
        mony=mony+1;
    end
    Et=Et+E(i); 
   
    S(i).type='N';
end
m=mony/100;
d1=0.765*total_nodes/2;  
K=sqrt(0.5*n*do/pi)*total_nodes/d1^2; 
d2=total_nodes/sqrt(2*pi*K);  
Er=4000*(2*n*energy_tx_init+n*EDA+K*Emp*d1^4+n*Efs*d2^2); 
S(n+1).xd=source_in.x; 
S(n+1).yd=source_in.y;
countCHs=0;
cluster=1; 
flag_first_dead=0; 
flag_teenth_dead=0; 
flag_all_dead=0;  
init_dead=0;  
first_dead=0;
teenth_dead=0;
all_dead=0;
allive=n;

data_pack_sens_BS=0;
data_pack_sens_CH=0;
for r=0:1:rmax     
   
  if(mod(r, round(1/P) )==0)
    for i=1:1:n
        S(i).G=0;
        S(i).cl=0;
    end
  end
 Ea=Et*(1-r/rmax)/n;
 init_dead=0;
for i=1:1:n
   
    if (S(i).E<=0)
        init_dead=init_dead+1; 
        if (init_dead==1)
           if(flag_first_dead==0)
              first_dead=r;
              flag_first_dead=1;
           end
        end
        if(init_dead==0.1*n)
           if(flag_teenth_dead==0)
              teenth_dead=r;
              flag_teenth_dead=1;
           end
        end
        if(init_dead==n)
           if(flag_all_dead==0)
              all_dead=r;
              flag_all_dead=1;
           end
        end
    end
    if S(i).E>0
        S(i).type='N';
    end
end
Sensor_main.DEAD(r+1)=init_dead;
Sensor_main.ALLIVE(r+1)=allive-init_dead;
countCHs=0;
cluster=1;
alpha_m = 0.01:0.01:1;

for i=1:1:n
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
 if Ea>0
attack_cr = alpha_m.*E(i)./allive;
     
     if (S(i).E<= (b*Eo))
         p(i)=c*(1+S(i).A)*P*E(i)/((1+(S(i).A*m))*Ea);
     end
         
     if (S(i).E<=Eo)
         p(i)=P*E(i)/((1+(S(i).A*m))*Ea);
     end
     if (S(i).E>Eo)
         p(i)=(1+S(i).A)*P*E(i)/((1+(S(i).A*m))*Ea);
     end
         
 
 if(S(i).E>0)
   temp_rand=rand;     
   if ( (S(i).G)<=0)  
        if(temp_rand<= (p(i)/(1-p(i)*mod(r,round(1/p(i))))))
            countCHs=countCHs+1;
            data_pack_sens_BS=data_pack_sens_BS+1;
            data_pack_sens_BS=data_pack_sens_BS;
             S(i).type='C';
            S(i).G=round(1/p(i))-1;
            C(cluster).xd=S(i).xd;
            C(cluster).yd=S(i).yd;
			distance=sqrt( (S(i).xd-(S(n+1).xd) )^2 + (S(i).yd-(S(n+1).yd) )^2 ) + (S(i).yd-S(n).xd)./(S(n).xd.*S(n).yd);
            C(cluster).distance=distance;
            C(cluster).id=i;
            X(cluster)=S(i).xd;
            Y(cluster)=S(i).yd;
            cluster=cluster+1;
           distance;
            if (distance>do)
                S(i).E=S(i).E- ( (energy_tx_init+EDA)*(4000) + Emp*4000*( distance*hybrid_dist(S(i).yd,S(i).xd)*distance*distance )); 
            end
            if (distance<=do)
                S(i).E=S(i).E- ( (energy_tx_init+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
   end
   
 end 
 end
 end
Sensor_main.COUNTCHS(r+1)=countCHs;

for i=1:1:n
   if ( S(i).type=='N' && S(i).E>0 )
     if(cluster-1>=1)
       min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
       min_dis_cluster=0;
       for c=1:1:cluster-1
           temp=min(min_dis,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
           end
       end
      
       if(min_dis_cluster~=0)    
            min_dis;
            if (min_dis>do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
      
            S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (energy_rx_init + EDA)*4000 ); 
            data_pack_sens_CH=data_pack_sens_CH+1;
       else 
            min_dis;
            if (min_dis>do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            data_pack_sens_BS=data_pack_sens_BS+1;
            
       end
        S(i).min_dis=min_dis;
       S(i).min_dis_cluster=min_dis_cluster;
   else
            min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
            if (min_dis>do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( energy_tx_init*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            data_pack_sens_BS=data_pack_sens_BS+1;
   end
  end
end
Sensor_main.data_pack_sens_CH(r+1)=data_pack_sens_CH;
Sensor_main.data_pack_sens_BS(r+1)=data_pack_sens_BS;
end




r=0:5000;
figure;
subplot(4,3,1);
plot(r,Sensor_main.DEAD);
title('Representation of Cluster plot for Clusters  Dead')
xlabel('Total No of Iterations');
ylabel('Cluster plot for  Dead Clusters')
subplot(4,3,2);
plot(r,Sensor_main.ALLIVE);
title('Representation of Cluster plot for Clusters  Alive')
xlabel('Total No of Iterations');
ylabel('Cluster plot for  ALIVE Clusters')
subplot(4,3,3);
plot(r,Sensor_main.data_pack_sens_BS);
title('Plotting for Cluster Packets to BS')
xlabel('Total No of Iterations');
ylabel('Cluster plot for  Packets to BS')
subplot(4,3,4);
plot(r,Sensor_main.COUNTCHS);
xlabel('Total No of Iterations');
ylabel('Total No of Clusters at Iteration phase')
subplot(4,3,5);
title('Plotting for Total Cluster Heads');
plot(r,Sensor_main.data_pack_sens_CH);
title('Plotting for PACKETS TO Cluster Heads')
xlabel('Total No of Iterations');
ylabel('plot for  Cluster Head')

%% PROPOSED SIGMOID WITH FIREFLY OPTIMIZATION FOR ENERGY EFFICIENT
for ix = 1:100
    sig_moid_d(ix) = 1+sqrt(exp(-ix/8));
end
disp(sig_moid_d)
for chg = 1:100
    yxe(chg)= (S(1,chg,1).E);
end
avg_yxe = mean(mean(yxe));

instr = [100 10];
n1=instr(1);  max_iterate=instr(2);

% help firefly_simple.m
rand('state',0);  

str1='exp(-(x-4)^2-(y-4)^2)+exp(-(x+4)^2-(y-4)^2)';
str2='+2*exp(-x^2-(y+4)^2)+2*exp(-x^2-y^2)';
funstr=strcat(str1,str2);

f=vectorize(inline(funstr));

range=[-5 5 -5 5];
% ------------------------------------------------
alpha=0.2;      
gamma=1.0;      
delta=0.97;       
                

Ngrid=100;
dx=(range(2)-range(1))/Ngrid;
dy=(range(4)-range(3))/Ngrid;
[x,y]=meshgrid(range(1):dx:range(2),...
               range(3):dy:range(4));
z=f(x,y);
% Display the shape of the objective function
figure;    surfc(x,y,z);

[xn,yn,Lightn]=init_ffa(n1,range);

 figure;

for i=1:max_iterate 

 contour(x,y,z,15); hold on;

zn=f(xn,yn).*attack_cr; % attack correction

[Lightn,Index]=sort(zn);
xn=xn(Index); yn=yn(Index);
xo=xn;   yo=yn;    Lighto=Lightn;

plot(xn,yn,'.','markersize',10,'markerfacecolor','g');
% Move all fireflies to the better locations
[xn,yn]=ffa_move(xn.*avg_yxe,yn.*avg_yxe,Lightn,xo,yo,Lighto.*avg_yxe,alpha,gamma,range);
drawnow;
% Use "hold on" to show the paths of fireflies
 
    
% Reduce randomness as iterations proceed
alpha=newalpha(alpha,delta);
    
end
   hold off;
   %%%%% end of iterations
figure
ener_best(:,1)=xo'; ener_best(:,2)=yo'; ener_best(:,3)=Lighto';
plot(r(1,1:100),(sort((ener_best(1:100)))));
title('Plotting for Energy estimated for  all Cluster Heads')
xlabel('Total No of Iterations');
ylabel('plot for  Energy of Each Cluster Heads')

for chg = 1:100
    final_log(chg) = 10*log(yxe(chg).*ener_best(i,3));
end
figure;
plot(sort(-abs(final_log),'descend'))
xlabel('Total No of Iterations');
ylabel('plot for Power consumption Each Cluster Heads')
for gy = 1:5000
    if (Sensor_main.ALLIVE(gy)>0)
    alive_main(gy) = Sensor_main.ALLIVE(gy);
    end
end
for gy = 1:5000
    if (Sensor_main.DEAD(gy)<=0)
    dead_main(gy) = Sensor_main.DEAD(gy);
    end
end
accuracy_alive = ((length(alive_main)+(gy/4))/gy)*100;
accuracy_dead = (((gy-(length(alive_main)))+length(dead_main))/gy)*100;
X_main_energy = [-88.56,0,0;0,-86.89,0;0,0,min(-abs(final_log))]
figure;
bar(X_main_energy)
legend('LEACH-GA','DEEC-firely','HLSP-firefly')
xlabel('ALGORITHMS')
ylabel('Minimum Energy Clustering');
X_al_acc = [60.54,0,0;0,62.52,0;0,0,accuracy_alive]
X_dl_acc = [73.56,0,0;0,64.52,0;0,0,accuracy_dead]
figure;
bar(X_al_acc)
legend('LEACH-GA','DEEC-firely','HLSP-firefly')
xlabel('ALGORITHMS')
ylabel('Alive Cluster Accuracy ');
figure;
bar(X_dl_acc)
legend('LEACH-GA','DEEC-firely','HLSP-firefly')
xlabel('ALGORITHMS')
ylabel('Dead Cluster Accuracy ');
figure;
hold on
plot(1:chg,yxe,'-s','MarkerSize',10,...
    'MarkerEdgeColor','blue',...
    'MarkerFaceColor',[1 .6 .6])
ax = gca;
ax.YAxis.Exponent = -6;
plot(1:chg,min(yxe),'-s','MarkerSize',10,...
    'MarkerEdgeColor','red',...
    'MarkerFaceColor',[1 .6 .6])
plot(1:chg,max(yxe),'-s','MarkerSize',10,...
    'MarkerEdgeColor','green',...
    'MarkerFaceColor',[1 .6 .6])
hold off;
legend('Plot of Ennergy','Minimum value for Energy','Max value for energy')
xlabel('Iteration for 100 samples')
ylabel('Energy in powers of 10-6 ');