% data - dataset, a row is a data point; K - the number of clusters
function MST_Split(data, K)

W = squareform(pdist(data));
% The first round to construct MST
Pairs1 = minimumSpanningTree(W);  
PP = Pairs1;
ddata = data;

[~, numOfVertices] = Prune(Pairs1);
idx = numOfVertices>0;
data = data(idx,:);
W = squareform(pdist(data));
Pairs1 = minimumSpanningTree(W); 

N = size(W,1);
k=fix(sqrt(N));

% The second round to construct MST
for ii = 1: 2
    
    if ii == 1
        Pairs_B = Pairs1;
    else
        Pairs_B = Pairs2;
    end;
    
    for i = 1:size(Pairs_B, 1)
        W(Pairs_B(i, 1), Pairs_B(i, 2)) = inf;
        W(Pairs_B(i, 2), Pairs_B(i, 1)) = inf;    
    end; 

    if ii == 1
        Pairs2 = minimumSpanningTree(W);  
    else
        Pairs3 = minimumSpanningTree(W);  
    end;
end;

% determine the centroids
G = zeros(N,2);
Pairs = [Pairs1;Pairs2;Pairs3];
W = squareform(pdist(data));

for i = 1: size(Pairs, 1)
    G(Pairs(i, 1), 1) = G(Pairs(i, 1), 1) + 1;
    G(Pairs(i, 1), 2) = G(Pairs(i, 1), 2) + W(Pairs(i, 1), Pairs(i, 2));
    
    G(Pairs(i, 2), 1) = G(Pairs(i, 2), 1) + 1;
    G(Pairs(i, 2), 2) = G(Pairs(i, 2), 2) + W(Pairs(i, 1), Pairs(i, 2));    
end;

DATA = zeros(N, N);

for i = 1: size(Pairs, 1)
    DATA(DATA(N,Pairs(i,1)) + 1, Pairs(i,1)) = Pairs(i,2);
    DATA(N,Pairs(i,1)) = DATA(N,Pairs(i,1)) + 1;
    
    DATA(DATA(N,Pairs(i,2)) + 1, Pairs(i,2)) = Pairs(i,1);
    DATA(N,Pairs(i,2)) = DATA(N,Pairs(i,2)) + 1;
end;

ii = 1;
C = zeros(k,1);
while ii <= k
    maxNum = 0;
    AvgLen = 0;
    idx = 0;
    for i = 1: N
        if G(i, 1) > maxNum
            maxNum = G(i,1);
            AvgLen = G(i,2) / G(i, 1);
            idx = i;
        elseif G(i, 1)~= 0 && G(i, 1) == maxNum && G(i, 2) / G(i, 1) < AvgLen
            maxNum = G(i,1);
            AvgLen = G(i,2) / G(i, 1);
            idx = i;
        end;    
    end;

    G(idx, :) = 0;
    for i = 1: DATA(N, idx)
        G(DATA(i, idx),1) = G(DATA(i, idx),1) - 1;
        G(DATA(i, idx),2) = G(DATA(i, idx),2) - W(idx, DATA(i, idx));
    end;
    
    C(ii) = idx;
    ii = ii + 1;
end;
Centroid = data(C,:);

% kmeans
IDX = kmeans(data, k, 'start', Centroid);
cluster_color = ['rbgmcr'];
shape = ['+oxv<>d']; 
[IDXSorted, IDX_index] =  sortrows(IDX);

% Find the all edges in a subset
x = data;
Pairs = Pairs1;

N = size(x,1);
DATA = zeros(N, N);
Pairs_IDX = zeros(N, N);
VisitedFlag = zeros(N);

for i = 1: size(Pairs, 1)
    DATA(DATA(N,Pairs(i,1)) + 1, Pairs(i,1)) = Pairs(i,2);
    Pairs_IDX(DATA(N,Pairs(i,1)) + 1, Pairs(i,1)) = i;
    DATA(N,Pairs(i,1)) = DATA(N,Pairs(i,1)) + 1;
    
    DATA(DATA(N,Pairs(i,2)) + 1, Pairs(i,2)) = Pairs(i,1);
    Pairs_IDX(DATA(N,Pairs(i,2)) + 1, Pairs(i,2)) = i;
    DATA(N,Pairs(i,2)) = DATA(N,Pairs(i,2)) + 1;
end;

Groups = zeros(N,k);
j= 1;
for i = 1: k
    while IDXSorted(j) == i    
        idx = IDX_index(j);
        j = j + 1;        
        for m = 1: DATA(N,idx)
            pairs_idx = Pairs_IDX(m,idx);
            if IDX(Pairs(pairs_idx,1)) == i && IDX(Pairs(pairs_idx,2)) == i && VisitedFlag(pairs_idx)==0
                Groups(Groups(N, i) + 1,i) = pairs_idx;
                Groups(N, i) = Groups(N, i) + 1;
                VisitedFlag(pairs_idx)=1;
            end;
        end;
        if j>size(IDXSorted, 1)
            break;
        end;
    end;
end;
    
% Find maximum substree
j = 1;
for i = 1: k
    m = 1;
    SubGroup = zeros(1, N);
    while j <= N && IDXSorted(j) == i
        SubGroup(1, IDX_index(j)) = m;
        j = j + 1;
        m = m + 1;
    end;
    
    Temp = zeros(m, m - 1);
    for p = 1: m - 1
        Temp(1, p) = IDX_index(j - m + p);
        Temp(m, p) = 1;
    end;
    
    for h = 1: Groups(N, i)
        pairs_index = Groups(h,i);
        pair1 = Pairs(pairs_index,1);
        pair2 = Pairs(pairs_index,2);
        
        sub_idx1 = SubGroup(1, pair1);
        sub_idx2 = SubGroup(1, pair2);
        for q = 1: Temp(m, sub_idx2)
            SubGroup(1, Temp(q, sub_idx2)) = sub_idx1;
            id = Temp(m, sub_idx1) + 1;            
            Temp(id, sub_idx1) = Temp(q, sub_idx2);
            Temp(m, sub_idx1) = Temp(m, sub_idx1) + 1;
        end;
        Temp(:,sub_idx2) = 0;                
    end;
    
    maxN = 0;
    max_idx = 0;
    for h = 1: m - 1
        if Temp(m, h) > maxN
            maxN = Temp(m, h);
            max_idx = h;
        end;
    end;
    
    for h = 1: m - 1
        if h == max_idx
            continue;
        end;
        for q = 1: Temp(m, h)
            IDX(Temp(q, h)) = 0;
        end;
    end;
end;

% redistribute 
Idx = find(IDX == 0);
nn = size(Idx, 1);
G = zeros(nn,1);%% groups of points with IDX == 0

gNo = 1;
for i = 1: N - 1
    if IDX(Pairs(i, 1)) == 0 && IDX(Pairs(i, 2)) == 0
        ii1 = find(Idx == Pairs(i, 1));
        ii2 = find(Idx == Pairs(i, 2));
        if G(ii1) == 0 && G(ii2) == 0
            G(ii1) = gNo;
            G(ii2) = gNo;
            gNo = gNo + 1;
        elseif G(ii1) == 0 && G(ii2) ~= 0
            G(ii1) = G(ii2);
        elseif G(ii1) ~= 0 && G(ii2) == 0
            G(ii2) = G(ii1);
        else
            idxx = G == G(ii1);
            G(idxx) = G(ii2);
        end;                         
    end;
end;

MinDist = zeros(gNo,2);
for i = 1: N - 1
    if IDX(Pairs(i, 1)) == 0 && IDX(Pairs(i, 2)) ~= 0
        ii = find(Idx == Pairs(i, 1));
        if G(ii) == 0
            continue;
        end;
        if W(Pairs(i, 1), Pairs(i, 2)) < MinDist(G(ii), 2) || MinDist(G(ii), 2) == 0
            MinDist(G(ii), 1) = Pairs(i, 2);
            MinDist(G(ii), 2) = W(Pairs(i, 1), Pairs(i, 2));
        end;
    end;
    if IDX(Pairs(i, 1)) ~= 0 && IDX(Pairs(i, 2)) == 0
        ii = find(Idx == Pairs(i, 2));
        if G(ii) == 0
            continue;
        end;
        if W(Pairs(i, 1), Pairs(i, 2)) < MinDist(G(ii), 2) || MinDist(G(ii), 2) == 0
            MinDist(G(ii), 1) = Pairs(i, 1);
            MinDist(G(ii), 2) = W(Pairs(i, 1), Pairs(i, 2));
        end;
    end;   
end;    

for i = 1: gNo
    idx = find(G == i);
    if size(idx, 1) == 0
        continue;
    end;
    IDX(Idx(idx)) = IDX(MinDist(i, 1));
end;

%%To deal with G == 0
idx = find(G == 0);
for i = 1: size(idx, 1)
    idxx = 0;
    w = 0;
    for j = 1: N -1  
        if Pairs(j, 1) == Idx(idx(i)) 
            if w == 0 || W(Pairs(j,1),Pairs(j,2))< w
                idxx = Pairs(j,2);
                w = W(Pairs(j,1),Pairs(j,2));
            end;
        end;
        if Pairs(j, 2) == Idx(idx(i)) 
            if w == 0 || W(Pairs(j,1),Pairs(j,2))< w
                idxx = Pairs(j,1);
                w = W(Pairs(j,1),Pairs(j,2));
            end;
        end;
    end;
    IDX(Idx(idx(i))) = IDX(idxx);
end;

%%% a patch
for i = 1: N
    closeIdx = 0;
    for j = 1: DATA(N,i)
        
        if closeIdx == 0 || W(i, DATA(j,i))<W(closeIdx,i)
            closeIdx = DATA(j,i);
        end;
    end;
    if IDX(i)~=IDX(closeIdx) && DATA(N,i)==2
        IDX(i) = IDX(closeIdx);
    end;
end;
                        
figure('Position',[400 200 400 400]);axis equal;
for i = 1:size(Pairs1, 1)
    d1 = [data(Pairs1(i,1),1);data(Pairs1(i,2),1)];
    d2 = [data(Pairs1(i,1),2);data(Pairs1(i,2),2)];
    path = line(d1, d2);
    set(path,'Color',[0,0,0],'LineWidth',0.5,'LineStyle','-','Clipping','off');hold on; 
end; 
data1=data';
for j=1:k
    id = find(IDX==j);
    plot(data1(1,id),data1(2,id),'o', 'MarkerEdgeColor', cluster_color(mod(j, 6)+1), 'MarkerFaceColor',[1,1,1],'MarkerSize',5); axis equal; hold on;
    %text(data1(1, id(1)), data1(2, id(1)), num2str(j),'FontSize',10);axis equal; hold on; 
end;

% Merge
figure('Position',[600 200 400 400]);
for i = k: -1: K + 1
    [IDX, num1, num2] = MergeGroup(W, IDX, i, N, DATA, Pairs1,Pairs2,Pairs3,Pairs_IDX,data);
    clf;
    for j=1: i - 1
        id = find(IDX==j);
        plot(data1(1,id),data1(2,id),shape(mod(j, 7)+1), 'MarkerEdgeColor', cluster_color(mod(j, 6)+1),'MarkerSize',5); axis equal; hold on;
        if j == num1
            text(data1(1, id(1)), data1(2, id(1)), num2str(j),'FontSize',18);axis equal; hold on; 
        end;
        text(data1(1, id(1)), data1(2, id(1)), num2str(j),'FontSize',10);axis equal; hold on; 
    end;    
end;

% merge hairs
IDX_Final = zeros(size(PP,1)+1,1);
Idx1 = find(numOfVertices > 0);
for i = 1: size(IDX,1)
    IDX_Final(Idx1(i)) = IDX(i);
end;

Idx2 = find(numOfVertices == 0);
for i = 1: size(Idx2,1)
    idxx = find(PP(:,1) == Idx2(i));
    if size(idxx,1) > 0
        idxx = PP(idxx, 2);
        idxxx = Idx1 == idxx;
        IDX_Final(Idx2(i)) = IDX(idxxx);
        continue;
    end;
    idxx = PP(:,2) == Idx2(i);
    idxx = PP(idxx, 1);
    idxxx = Idx1 == idxx;
    IDX_Final(Idx2(i)) = IDX(idxxx);
end;

figure('Position',[600 200 400 400]);
data1 = ddata';
for j=1: K
    id = find(IDX_Final==j);
    plot(data1(1,id),data1(2,id),shape(mod(j, 7)+1), 'MarkerEdgeColor', cluster_color(mod(j, 6)+1),'MarkerSize',5); axis equal; hold on; 
end;    

%MergeGroup
function [IDX, num1, num2] = MergeGroup(W, IDX, k, N, DATA, Pairs1,Pairs2,Pairs3,Pairs_IDX,data)
% Find the group pairs to be considered for merging
Merge = zeros(k + 1,k);
[IDXSorted, IDX_index] =  sortrows(IDX);
j = 0;
s = 1;
for i = 1: k
    j = j + 1;
    while s <= N && IDXSorted(s) == j
        idx = IDX_index(s);
        for h = 1: DATA(N,idx)
            pair_idx = Pairs_IDX(h, idx);
            idx1 = Pairs1(pair_idx, 1);
            idx2 = Pairs1(pair_idx, 2);

            if IDX(idx1) ~= i
                Merge(Merge(k+1, i) + 1, i) = IDX(idx1);
                Merge(k+1, i) = Merge(k+1, i) + 1;
            elseif IDX(idx2) ~= i
                Merge(Merge(k+1, i) + 1, i) = IDX(idx2);
                Merge(k+1, i) = Merge(k+1, i) + 1;
            end;
        end;
        s = s + 1;
    end;       
end;

% Find the number of straddle edges and the sum of weights
NumEdge = zeros(k,k);
SumWeight = zeros(k,k);
ConnectedNumEdge = zeros(k,k);
ConnectedSumWeight = zeros(k,k);

densityN = fix(N/k *3);
DensitySpanIdx = zeros(k,k,densityN);
DensitySpan = zeros(k,k);

Pairs = [Pairs1;Pairs2;Pairs3];
Pairs_IDX = zeros(N,N);
for i = 1: size(Pairs, 1)
    Pairs_IDX(Pairs_IDX(N,Pairs(i,1)) + 1, Pairs(i,1)) = i;
    Pairs_IDX(N,Pairs(i,1)) = Pairs_IDX(N,Pairs(i,1)) + 1;
    
    Pairs_IDX(Pairs_IDX(N,Pairs(i,2)) + 1, Pairs(i,2)) = i;
    Pairs_IDX(N,Pairs(i,2)) = Pairs_IDX(N,Pairs(i,2)) + 1;
end;

for i = 1: size(Pairs, 1)
    idx1 = Pairs(i, 1);        
    idx2 = Pairs(i, 2);
    
    if IDX(idx1) ~= IDX(idx2)
        %straddling
        NumEdge(IDX(idx1), IDX(idx2)) = NumEdge(IDX(idx1), IDX(idx2)) + 1;
        NumEdge(IDX(idx2), IDX(idx1)) = NumEdge(IDX(idx2), IDX(idx1)) + 1;
        
        SumWeight(IDX(idx1), IDX(idx2)) = SumWeight(IDX(idx1), IDX(idx2)) + W(idx1,idx2);
        SumWeight(IDX(idx2), IDX(idx1)) = SumWeight(IDX(idx2), IDX(idx1)) + W(idx1,idx2);
        
        ii = DensitySpanIdx(IDX(idx1), IDX(idx2), densityN) + 1;
        DensitySpanIdx(IDX(idx1), IDX(idx2), densityN) = ii;
        DensitySpanIdx(IDX(idx1), IDX(idx2), ii) = idx1;

        ii = DensitySpanIdx(IDX(idx2), IDX(idx1), densityN) + 1;
        DensitySpanIdx(IDX(idx2), IDX(idx1), densityN) = ii;
        DensitySpanIdx(IDX(idx2), IDX(idx1), ii) = idx2;          
        
        %connected
        for j = 1: Pairs_IDX(N, idx1)
            pairs_idx = Pairs_IDX(j, idx1);
            if IDX(Pairs(pairs_idx, 1)) == IDX(Pairs(pairs_idx, 2))
                ConnectedNumEdge(IDX(idx1),IDX(idx2)) = ConnectedNumEdge(IDX(idx1),IDX(idx2)) + 1;
                ConnectedSumWeight(IDX(idx1),IDX(idx2)) = ConnectedSumWeight(IDX(idx1),IDX(idx2)) + W(Pairs(pairs_idx, 1),Pairs(pairs_idx, 2));              
            end;
        end;
        for j = 1: Pairs_IDX(N, idx2)
            pairs_idx = Pairs_IDX(j, idx2);
            if IDX(Pairs(pairs_idx, 1)) == IDX(Pairs(pairs_idx, 2))
                ConnectedNumEdge(IDX(idx2),IDX(idx1)) = ConnectedNumEdge(IDX(idx2),IDX(idx1)) + 1;
                ConnectedSumWeight(IDX(idx2),IDX(idx1)) = ConnectedSumWeight(IDX(idx2),IDX(idx1)) + W(Pairs(pairs_idx, 1),Pairs(pairs_idx, 2));               
            end;
        end;
    end;
end;

for i = 1: k
    for j = 1: k
        Idx = DensitySpanIdx(i,j,1:DensitySpanIdx(i,j,densityN));
        if size(Idx,3) == 0
            continue;
        end;
        
        maxL = max(pdist(data(Idx,:)));
        if size(maxL,2) == 0 
            continue;
        end;        
        DensitySpan(i,j) = maxL;
    end;
end;

% Compute the attractors
Attractor = zeros(k,k);

% intracluster similarity;
IntraSimi = zeros(k,2);
for i = 1: k
    IDX_S = find(IDX == i);
    [no,weight] = bisect(IDX_S, data,k);    
    IntraSimi(i, 1) = no;
    IntraSimi(i, 2) = weight;
end;

for i = 1: k
    size_i = size(find(IDX == i),1);
    for j = 1: Merge(k+1, i)
        num = NumEdge(i,Merge(j,i));
        avgWeightStraddle = SumWeight(i,Merge(j,i))/NumEdge(i,Merge(j,i));
        if ConnectedNumEdge(i,Merge(j,i)) == 0
            avgWeightConnected1 = 1;
        else
            avgWeightConnected1 = ConnectedSumWeight(i,Merge(j,i))/ConnectedNumEdge(i,Merge(j,i));
        end;
        
        if ConnectedNumEdge(Merge(j,i),i) == 0
            avgWeightConnected2 = 1;
        else
            avgWeightConnected2 = ConnectedSumWeight(Merge(j,i),i)/ConnectedNumEdge(Merge(j,i),i);
        end;        
        
        avgWeightConnected = min(avgWeightConnected1,avgWeightConnected2);
        density_span = max(DensitySpan(i,Merge(j,i)), DensitySpan(Merge(j,i),i));
        
        size_j = size(find(IDX == j),1);
        att = num / min(size_i,size_j); 
        att = att * (1/(avgWeightStraddle/avgWeightConnected)) * density_span;    
        
        no1 = IntraSimi(i,1);
        weight1 = IntraSimi(i,2);
        no2 = IntraSimi(Merge(j,i),1);
        weight2 = IntraSimi(Merge(j,i),2);
                
        factor = no1*no2/(no1+no2);
        factor = factor * (no1+no2);
        factor = (sqrt(weight1/no1 * weight2/no2)) /(weight1/no1 + weight2/no2)/factor;

        no3 = max(no1, no2);
        factor1 = min(num, no3)/max(num, no3);
        
        att = att * factor;
        %att = att * factor1;   

        Attractor(i,Merge(j,i)) = att;
        Attractor(Merge(j,i),i) = att;
    end;
end;

% Merge and Update 
maxAtt = 0;
idx1 = 0;
idx2 = 0;
for i = 1: k
    for j = 1: k
        if Attractor(i,j) > maxAtt
            maxAtt = Attractor(i,j);
            idx1 = i;
            idx2 = j;
        end;
    end;
end;

num1 = idx1;
num2 = idx2;

% update IDX
idxMax = 0;
idxMin = 0;
if idx1 > idx2
    idxMax = idx1;
    idxMin = idx2;   
else
    idxMax = idx2;
    idxMin = idx1;
end;

for i = 1: N
    if IDX(i) == idxMax
        IDX(i) = idxMin;
    end;
end;

for i = 1: N
    if IDX(i) > idxMax
        IDX(i) = IDX(i) - 1;
    end;
end;


% Kruskal's algorithm 
function Pairs = minimumSpanningTree(W) 
% the third column of path is the representative element flag
Pairs = zeros(size(W, 1) - 1, 3);
D = zeros(size(W, 1) * (size(W, 1)  - 1) / 2, 1);
I = zeros(size(D,1), 2);
m = 1;
n = 1;

for i = 1: size(D,1)

    I(i, 1) = m + 1;
    I(i, 2) = n;
    D(i) = W(m+ 1, n);
    
    if n == m
        m = m + 1;
        n = 1;
    else
        n = n + 1;
    end;
end;

D = [D, I];
D = sortrows(D);

numOfPairsIn = 0;
for i = 1: size(D, 1)
    rep1 = repElement(Pairs, D(i, 2));
    rep2 = repElement(Pairs, D(i, 3));
    
    if rep1 ~= rep2 || (rep1 < 0 && rep2 < 0)
        numOfPairsIn = numOfPairsIn + 1;
        Pairs(numOfPairsIn, 1) = D(i, 2);
        Pairs(numOfPairsIn, 2) = D(i, 3);
        
        if rep1 < 0 && rep2 <0
            Pairs(numOfPairsIn, 3) = D(i, 2);
        elseif rep1 < 0 && rep2 > 0
            Pairs(numOfPairsIn, 3) = rep2;
        elseif rep1 > 0 && rep2 < 0
            Pairs(numOfPairsIn, 3) = rep1;
        else                                % if both rep1 and rep2 are greater than 0
            Pairs(numOfPairsIn, 3) = rep1;
            for j = 1: size(Pairs, 1)       % combine two subtrees 
                if Pairs(j,3) == rep2
                    Pairs(j,3) = rep1;
                end;
            end;
        end;
    end;
    
    if numOfPairsIn >= size(W, 1) - 1
        Pairs(:,3) = 0;
        break;
    end;
end;

Pairs = Pairs(:,1:2);

% determine the representative element
function rep = repElement(Pairs, node)
rep = -1;
for i = 1: size(Pairs, 1)
    if Pairs(i, 1) == node || Pairs(i, 2) == node
        rep = Pairs(i, 3);
        return;
    end;
end;

% Pruning hairs
% cut the hairs,until two iterations have cut off the same number of hairs
function [Pairs, numOfVertices] = Prune(Pairs)
n = size(Pairs, 1);
numOfVertices = zeros(n+1,1);
numChanged = 0;

iii = 1;
while true
    numChanged1 = 0;
    for i = 1: n
        if Pairs(i,1) == 0
            continue;
        end;
        numOfVertices(Pairs(i,1)) = numOfVertices(Pairs(i,1)) + 1;
        numOfVertices(Pairs(i,2)) = numOfVertices(Pairs(i,2)) + 1;
    end;

    for i = 1: n
        if Pairs(i,1) == 0 
            continue;
        end;
        if numOfVertices(Pairs(i,1)) == 1 || numOfVertices(Pairs(i,2)) == 1
            Pairs(i,:) = 0;
            numChanged1 = numChanged1 + 1;
        end;
    end;
    
    iii = iii + 1;
    if iii == 2 %numChanged == numChanged1
        break;
    end;
    numChanged = numChanged1;
    numChanged1 = 0;
    numOfVertices(:) = 0;
end;

numOfVertices(:) = 0;
Pairs = Pairs(Pairs(:,1)>0,:);
for i = 1: size(Pairs,1)
    numOfVertices(Pairs(i,1)) = numOfVertices(Pairs(i,1)) + 1;
    numOfVertices(Pairs(i,2)) = numOfVertices(Pairs(i,2)) + 1;
end;

% bisect
function [no, weight] = bisect(IDX_S, data,k)

D = data(IDX_S,:);
W = squareform(pdist(D));
W1 = W;
Pairs1 = minimumSpanningTree(W); 

for i = 1: size(Pairs1, 1)
    W(Pairs1(i,1),Pairs1(i,2)) = inf;
    W(Pairs1(i,2),Pairs1(i,1)) = inf;
end;
Pairs2 = minimumSpanningTree(W); 

for i = 1: size(Pairs2, 1)
    W(Pairs2(i,1),Pairs2(i,2)) = inf;
    W(Pairs2(i,2),Pairs2(i,1)) = inf;
end;
Pairs3 = minimumSpanningTree(W); 

nn = size(Pairs1,1);
G = zeros(nn + 1, 1);
idx = 0;
bRation = 0;

for i = 1: nn
    G_B = zeros(nn + 1, 1);
    no = 1;
    for j = 1: nn
        if j == i 
            continue;
        end;
        
        if G_B(Pairs1(j,1)) == 0 && G_B(Pairs1(j,2)) == 0
            G_B(Pairs1(j,1)) = no;
            G_B(Pairs1(j,2)) = no;
            no = no + 1;
        elseif G_B(Pairs1(j,1)) ~= 0 && G_B(Pairs1(j,2)) == 0
            G_B(Pairs1(j,2)) = G_B(Pairs1(j,1));
        elseif G_B(Pairs1(j,1)) == 0 && G_B(Pairs1(j,2)) ~= 0
            G_B(Pairs1(j,1)) = G_B(Pairs1(j,2));
        elseif G_B(Pairs1(j,1)) ~= 0 && G_B(Pairs1(j,2)) ~= 0
            idxx1 = Pairs1(j,1);
            idxx2 = Pairs1(j,2);
            if G_B(idxx1) < G_B(idxx2)
                G_B(G_B == G_B(idxx2)) = G_B(idxx1);
            else
                G_B(G_B == G_B(idxx1)) = G_B(idxx2);
            end;
        end;        
    end;
    no1 = size(find(G_B == 1),1);
    no1 = min(no1, nn -no1 + 1)/max(no1, nn -no1 + 1);
    if bRation == 0 || no1 > bRation
        G = G_B;
        idx = i;
        bRation = no1;
    end;
end;

no = 1;
weight = W1(Pairs1(idx,1),Pairs1(idx,2));
for i = 1: nn
    if G(Pairs2(i, 1)) ~= G(Pairs2(i,2))
        no = no + 1;
        weight = weight + W1(Pairs2(i,1), Pairs2(i,2));
    end;
    if G(Pairs3(i, 1)) ~= G(Pairs3(i,2))
        no = no + 1;
        weight = weight + W1(Pairs3(i,1), Pairs3(i,2));
    end;
end;
