SystemEvolution_energy(labels,Dist,dc,N,Ns,verbose)

function [inter,inter2,intra,intra2,N,Ns,nCluster] = SystemEvolution_energy(labels,Dist,dc,N,Ns,verbose)
% Energy computations for System Evolution to Estimate the Number of Clusters


%(1) calculating partition & merge energies
if 1 
inter=zeros(N,4);     inter2=zeros(N,4);
intra=zeros(N,4);     intra2=zeros(N,4); 
inter3=zeros(N,4);   intra3=zeros(N,4);
nCluster=zeros(1,N); 
end
plots=0;                    % plotting intermediate solutions ?

for i = Ns
   if verbose
   fprintf('\n = = system Evolution runs at k= %d ', i);
   end
   Q = ind2cluster(labels(:,i));               % fetching labels of clusters
   [g,h,k,ters] = close_cluster(Dist,Q,dc,1); % find closest-cluster candidates
   if k == 0                    % if less than 4 points in a cluster
     N = i-1;
     k = Ns < i;
     Ns = Ns(k);
     if verbose
     fprintf('\n # # Stop at k= %d ! a cluster has points < 4', i);
     end
     break; 
   end
   k = length(g);
   Ep=[];
   Em = [];
   
   for j = 1: k
     % specifying border and central regions of each cluster
     [A1,A2,B1,B2,A3,B3] = cluster_model(Dist,Q{g(j)},Q{h(j)});
     % calculating energies between each pair of candidate closest clusters
     Ep(j,:) = energy_compute(Dist,A1,A2,B1,B2,A3,B3,dc);
     if plots
       clf; plotdata_bylabels(data,labels(:,i),pc,subp,'nb');
       plotdata_bylabels(data,labels(:,i),5,subp,'co',A1,A2,B1,B2,A3,B3);
     end
   end
   [ind,g1] = sort(Ep(:,4)); 
   g2 = g1(1);     % minimal energy difference among k pairs of candidates
   inter(i,:) = Ep(g2,:);
   nCluster(i)=min([length(Q{g(g2)}) length(Q{h(g2)})]);
   if i > 2
      inter2(i,:) = Ep(g1(2),:);
      inter3(i,:) = Ep(g1(3),:);
   end
   
   if inter(i,4) < 0 && (inter2(i,4) < 0 || i < 4) % overlapping case
      [g,h] = close_cluster(Dist,Q,dc,2); % find closest-cluster candidates
      for j = 1: k
         [A1,A2,B1,B2,A3,B3] = cluster_model(Dist,Q{g(j)},Q{h(j)});
         [Ep(j,:),Em(j,:)] = energy_compute(Dist,A1,A2,B1,B2,A3,B3,dc+3);
         if plots
         clf; plotdata_bylabels(data,labels(:,i),5,subp,'co',A1,A2,B1,B2,A3,B3);
         end
      end
      ind = find(Ep(:,4)<0);
      [g1, km] = sort(Em(ind,4));
      km = ind(km);
      ind2 = find(Ep(:,4)>0);
      [g1, kn] = sort(Em(ind2,4));
      kn = ind2(kn);
      g1 = [km; kn];
      % finding those with Ep<0 but Em>0 from k pairs of candidates
      ind2 =[];
      for  j = 1: k 
        if Em(g1(j),4) < 0 && Ep(g1(j),4) > 0
           ind2 = [ind2 j];
        end
      end
      g1(ind2) = [];
      intra(i,:) = Em(g1(1),:);
      ind = length(g1);
      if i > 2 && ind > 1 && min([length(Q{g(g1(2))}) length(Q{h(g1(2))})]) >10
         intra2(i,:) = Em(g1(2),:);
         if ind > 2 && min([length(Q{g(g1(3))}) length(Q{h(g1(3))})]) >10
            intra3(i,:) = Em(g1(3),:);
         end
      end
   elseif inter(i,4) < 0 % if second smallest distance inter2(i,4) > 0
      ind = [];
      if i > 3
         ind2 = [];
         Ep = ters(h(g2),:);
         [Ep, km] = sort(Ep);
         ind(:,2) = km(2:i)';
         ind(:,1) = h(g2)';
         kn = find(ind(:,2)==g(g2));
         Ep = ters(g(g2),:);
         Ep(h(g2)) = 0;
         [Ep, km] = sort(Ep);
         ind2(:,2) = km(3:i)';
         ind2(:,1) = g(g2)';
         ind = [ind;ind2]; % checking all the clusters from twin-clusters
      else
         ind = [g' h'];
         kn = 3;
      end
      Ep = [];
      Em = [];
      for j = 1:size(ind,1)
         [A1,A2,B1,B2,A3,B3]=cluster_model(Dist,Q{ind(j,1)},Q{ind(j,2)});
         [Ep(j,:),Em(j,:)] = energy_compute(Dist,A1,A2,B1,B2,A3,B3,dc+3);
         if plots
         clf; plotdata_bylabels(data,labels(:,i),5,subp,'co',A1,A2,B1,B2,A3,B3);
         end
         if Ep(j,4) <= 0 && Em(j,4) <= 0 % if any two clusters are inseparable
            intra(i,:) = Em(j,:);
            g2 = 0;
            break;
         end
      end
      if g2 > 0
         Ep = [Ep; Ep(kn,:)];
         Em = [Em; Em(kn,:)];
         ind2 = Ep(:,4);
         ind2 = find(ind2>0);
         ters = Em(:,4);
         ters(ind2) = 2;
         [ko g1] = min(ters(1:i-1));
         [km g2] = min(ters(i:i+i-2));
         g2 = g2+i-1;
         if ko <= km
            intra(i,:) = Em(g1,:);
         else
            intra(i,:) = Em(g2,:);
         end
      end
   end
end