[inter,intra]=energy_pearson(Dmatrix,A1,A2,B1,B2,A3,B3,chois,fine)

function [inter,intra]=energy_pearson(Dmatrix,A1,A2,B1,B2,A3,B3,chois,fine)
% computing energies (average walking distance) within each group A1,B1,...

k = mod(chois,3); 
kmin = 6;
% computing energies for well-separated twin-clusters
Dm = Dmatrix;
ch = 1;                        % using max(dista distb)
Oa =[]; 
Ob= [];
 
% finding average walk distance within every A1,B1,...
nA = length(A1);
A1o = A1;
if nA > 10
   [A1o,Oa,Dm] = remove_outlier(Dmatrix,A1,nA);
end
distA1 = walk_distance(Dmatrix,A1o,k,ch);
Ra = remove_outlier(Dmatrix,A2,length(A2));
distA2 = walk_distance(Dmatrix,Ra,k,ch);
distA12 = walk_distance(Dmatrix,[A1;A2],k,ch);
if length(A3) > 0
   A3 = remove_outlier(Dmatrix,A3,length(A3));
   distA3 = walk_distance(Dmatrix,A3,k,ch);
   distA23 = walk_distance(Dmatrix,[A2;A3],k,ch);
   distA2h = max([distA2 2*(distA3*distA23)/(distA3+distA23)]); 
   distA2l = min([distA2 distA3 distA23]); 
else
   distA2h = distA2; 
   distA2l = distA2; 
   if length(A2) < kmin
      distA2l = max([distA12 distA2]); 
   end
end

nB = length(B1);
B1o = B1;
if nB > 10
   [B1o,Ob,Dm] = remove_outlier(Dmatrix,B1,nB);
end
distB1 = walk_distance(Dmatrix,B1o,k,ch);
Rb = remove_outlier(Dmatrix,B2,length(B2));
distB2 = walk_distance(Dmatrix,Rb,k,ch);
distB12 = walk_distance(Dmatrix,[B1;B2],k,ch);
if length(B3) > 0
   B3 = remove_outlier(Dmatrix,B3,length(B3));
   distB3 = walk_distance(Dmatrix,B3,k,ch);
   distB23 = walk_distance(Dmatrix,[B2;B3],k,ch);
   distB2h = max([distB2 2*(distB3*distB23)/(distB3+distB23)]); 
   distB2l = min([distB2 distB3 distB23]); 
else
   distB2h = distB2; 
   distB2l = distB2;
   if length(B2) < kmin 
      distB2l = max([distB2 distB12]);
   end
end
  
% finding nearest elements between A1 & B1
na = 2*fix(sqrt(nA));
Ra = border_points(Dmatrix,A1,B1,na);
nb = 2*fix(sqrt(nB));
Rb = border_points(Dmatrix,B1,A1,nb);
Ra = setdiff(Ra,Oa);
Rb = setdiff(Rb,Ob);
if na > 5 && nb > 5
   [Ra,Rb] = remove_closest(Dmatrix,Ra,Rb);
end
na = length(Ra);
nb = length(Rb);

% computing border distance between A1 & B1
if na == 1
  distab = [Dmatrix(Ra,Rb) Dmatrix(Ra,Rb)];
elseif nb == 1
  distab = [Dmatrix(Ra,Rb)' Dmatrix(Ra,Rb)'];
else
  distab = [min(Dmatrix(Ra,Rb)) min(Dmatrix(Rb,Ra))];
end
[distab,cut] = sort(-distab);
cut = max([fix(sqrt(na+nb))-1 1]);
distab = -distab;
nb = min([na+nb+1-cut na+nb]);
cut = distab(cut:nb);
distAB = mean(cut);
distAB1 = walk_distance(Dmatrix,[Ra;Rb],k,10); %10: keep bigest route
cut = distAB/distAB1;
if  cut > 1 && cut < 3   % if twin-clusters are close [A1o;B1o]
   distAB = 0.5*(distAB+distAB1);
end

% output
distA = max([distA1 distA2h distA12]); %distA2
distB = max([distB1 distB2h distB12]);
cut = max([distA distB]);
cut = (distAB-cut);
inter = [distA distB distAB cut];
if chois < 3  %|| cut > 0.01 for a look of computations in the follwoing
   intra = inter;
   return;
end

% computing energies for overlapping twin-clusters
ch = 2;                        % using mean(dista distb)
na = round(0.5*nA);
Ra = border_points(Dmatrix,A1,B1,na);
nb = round(0.5*nB);
Rb = border_points(Dmatrix,B1,A1,nb);
distA1b = walk_distance(Dmatrix,[A1;Rb],k,ch); % across A1-B1
distB1a = walk_distance(Dmatrix,[B1;Ra],k,ch);

distAi = 2*(distA1*distA1b)/(distA1+distA1b);
distBi = 2*(distB1*distB1a)/(distB1+distB1a);
distAj = 2*(distAi*distA2l)/(distAi+distA2l); % integrate A1-B2
distBj = 2*(distBi*distB2l)/(distBi+distB2l);
if fine
   distA = min([distAj distA12]);
   distB = min([distBj distB12]);
   if length(A2) < kmin
      distA = max([distAj distA12]);
   end
   if length(B2) < kmin
      distB = max([distBj distB12]);
   end
else
   distA = max([distAj distA12]);
   distB = max([distBj distB12]);
end


% calculating cross region F (A1-B1)
distAB1 = walk_distance(Dmatrix,[A1o;B1o],k,ch);
[Ra,Rb] = regionF_find(Dmatrix,A1o,B1o,distA1,distB1,0.9);
distAB2 = walk_distance(Dmatrix,[Ra;Rb],k,ch);
[Ra,Rb] = regionF_find(Dmatrix,A1o,B1o,distA1,distB1,0.8);
distAB3 = walk_distance(Dmatrix,[Ra;Rb],k,ch);
[Ra,Rb] = regionF_find(Dmatrix,A1o,B1o,distA1,distB1,0.7);
distAB4 = walk_distance(Dmatrix,[Ra;Rb],k,ch);

na = sort([distAB1 distAB2 distAB3 distAB4]);
na = mean(na(3:4)); 
nb = max([distAB1 distAB2 distAB3]);
cut = max([distA distB]);
Ra = na-cut; %(na-cut)/cut
Rb = nb-cut; %(nb-cut)/cut

% output
if  na > nb
intra = [distA distB na Ra];
else
intra = [distA distB nb Rb];
end