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

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

k = mod(chois,3); 
% computing energies for well-separated twin-clusters
ch = 1;                        % using max(dista distb)
Oa =[]; 
Ob= [];
nA = length(A1);
nB = length(B1);
A1o = A1;
B1o = B1;
if nA > 10
   [A1o, Oa, Dmatrix] = remove_outlier(Dmatrix,A1,nA);
end
if nB > 10
   [B1o, Ob, Dmatrix] = remove_outlier(Dmatrix,B1,nB);
end
Dm = Dmatrix;

% finding average walk distance within every A1,B1,...
distA1 = walk_distance(Dmatrix,A1,k,ch);
distA2 = walk_distance(Dmatrix,A2,k,ch);
distA12 = walk_distance(Dmatrix,[A1;A2],k,ch);
if length(A3) > 0
   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; 
end

distB1 = walk_distance(Dmatrix,B1,k,ch);
distB2 = walk_distance(Dmatrix,B2,k,ch);
distB12 = walk_distance(Dmatrix,[B1;B2],k,ch);
if length(B3) > 0
   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; 
end
  
% finding nearest elements between A1 & B1
na = 2*fix(sqrt(nA));
Ra = border_points(Dm,A1,B1,na);
nb = 2*fix(sqrt(nB));
Rb = border_points(Dm,B1,A1,nb);
Ra = setdiff(Ra,Oa);
Rb = setdiff(Rb,Ob);
if na > 5 && nb > 5
   [Ra,Rb] = remove_closest(Dm,Ra,Rb);
end
na = length(Ra);
nb = length(Rb);

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

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

% computing energies for overlapping twin-clusters
ch = 2;                         % using mean(dista distb)
na = round(0.5*nA);
Ra = border_points(Dm,A1,B1,na);
nb = round(0.5*nB);
Rb = border_points(Dm,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);
distA = 2*(distAi*distA2l)/(distAi+distA2l); % integrate A1-B2
distB = 2*(distBi*distB2l)/(distBi+distB2l);
distA = max([distA distA2 distA2h]); 
distB = max([distB distB2 distB2h]);

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

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

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