[A1,A2,B1,B2,A3,B3]=cluster_model(Dmatrix,Ua,Ub)

function [A1,A2,B1,B2,A3,B3]=cluster_model(Dmatrix,Ua,Ub)
% Two-cluster model and twin-clusters
% input: Ua/Ub - labels for cluster A & B, Dmatrix: distance/similarity matrix
% output: border region A1 & B1 are close each other; A2/B2 closer to A1/B1

% nA/nB be total number of elements of cluster A/B
nA = length(Ua); 
nB = length(Ub); 

% each cluster is cut into m parts, then nA1=nA/m, nB1=nB/m
cutA = 3;
cutB = 3;

if nA > 49
   cutA = fix(0.6*sqrt(nA));
end
if nB > 49
   cutB = fix(0.6*sqrt(nB));
end

% m=2 is recommended for nA<21
if nA < 21
   mA = fix(0.5*nA); 
else
   mA = fix(nA/cutA)+1; 
end
if nB < 21
   mB = fix(0.5*nB); 
else
   mB = fix(nB/cutB)+1; 
end
if mA > 11
   cut = round(0.1*mA);
   mA = mA -cut;
end
if mB > 11
   cut = round(0.1*mB);
   mB = mB -cut;
end

% finding 2 border regions A1/A2 and B1/B2
A1 = border_points(Dmatrix,Ua,Ub,mA);
A = setdiff(Ua,A1);
A2 = border_points(Dmatrix,A,A1,mA);
B1 = border_points(Dmatrix,Ub,Ua,mB);
B = setdiff(Ub,B1);
B2 = border_points(Dmatrix,B,B1,mB);

% finding one more region if cut parts >3
A3 = [];
B3 = [];
if cutA > 3
  A = setdiff(A,A2);
  A3 = border_points(Dmatrix,A,A2,mA);
end

if cutB > 3
  B = setdiff(B,B2);
  B3 = border_points(Dmatrix,B,B2,mB);
end
%[A1,B1]=findnear(Dmatrix,Ua,mA,Ub,mB,2);