function [Nclust, LabelClust, assig] = clustering_lzl(dis_matrix,percent)
%%
% Input  - dis_matrix : distance matrix in which lower value means high similarity
%          percent    : a parameter that determines the density of data point (0 - 1)
%          out        : output folder
% Output - Nclust     : number of clusters
%          LabelClust £ºcluster centers
%          assig      £ºassigment of each data point 

% reference:  Clustering by fast search and find of density peaks, Science, 2014
% 2016/11/28, UESTC, LZL
%%
N = size(dis_matrix,1); % 
dis_matrix(1:N+1:end) = 0; % 

%% find threshold of dc
% method 1.
Minv = min(dis_matrix(:));
Maxv = max(dis_matrix(:));
for i=Minv : (Maxv-Minv)/1000 : Maxv
    A = dis_matrix < i & dis_matrix > 0;
    A = mean(sum(A,2));
    if A >= ceil(N*percent)
        break;
    end
end
dc = i;
clear i A Minv Maxv

% % method 2. 
% Ind = find(triu(ones(N),1));
% xx  = dis_matrix(Ind);
% xx  = sort(xx);
% dc = xx(ceil(length(Ind)*percent));
% clear Ind xx

%% find RUO £¨density of each data point£©
RUO = zeros(N,1);
for i=1:N
    A = dis_matrix(i,setdiff(1:N,i));
    A(A == 0) = [];
    for iInd = 1:length(A)
        RUO(i) = RUO(i) + exp(-(A(iInd)/dc)*(A(iInd)/dc)); % gaussion kernal
    end
    clear A iInd
end
clear i

%% find delta £¨the smallest distance between each data point and those points with high density values£©
DELTA = zeros(N,1);
[~,RUOind] = sort(RUO,'descend');
for i=1:length(RUOind)
    if i==1
        DELTA(RUOind(i)) = max(dis_matrix(:));
    else
        A = dis_matrix(RUOind(i),RUOind(1:i-1));
        DELTA(RUOind(i)) = min(A);
        clear A
    end
end
clear RUOind i 
% 

%% figuring
scrsz = get(0,'ScreenSize');
H = figure('Position',[200 100 scrsz(4)*1.4 scrsz(4)/2]);

subplot(1,2,1);
plot(RUO,DELTA,'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k'); hold on;
title ('Decision Graph','FontSize',10.0);
xlabel ('\rho');
ylabel ('\delta');

subplot(1,2,2);
% figure;
gamma = RUO.*DELTA;
[A, Aind] = sort(gamma,'descend');
xrange = 1:N;
plot(xrange,A,'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k'); hold on;
title (['\rho',' x ','\delta'],'FontSize',10.0);
xlabel('n');
ylabel('\gamma');
clear A

Nclust = inputdlg({'Number of clusters'},'Parameters',1,{'1'});
Nclust = str2double(Nclust);
fprintf('The number of clusters are £º %i \n',Nclust);

LabelClust = zeros(Nclust,1); % 
cmap       = colormap;
cmap_use   = [];
RUO_temp   = RUO;
DELTA_temp = DELTA;
gamma_temp = gamma;
for i = 1:Nclust
    subplot(1,2,1)
    rect     = getrect(H);
    rhomin   = rect(1);
    deltamin = rect(2);
    
    D        = sqrt((RUO_temp - rhomin).^2 + (DELTA_temp - deltamin).^2);
    [~,minInd] = min(D);
    LabelClust(i) = minInd;
    ic       = int8((i*64.)/(Nclust*1.));
    
    subplot(1,2,1)
    hold on;
    plot(RUO_temp(minInd),DELTA_temp(minInd),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
    
    subplot(1,2,2)
    hold on;
    plot(xrange(Aind == minInd),gamma_temp(minInd),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
    
    cmap_use = [cmap_use;cmap(ic,:)];
    
    RUO_temp(minInd) = Inf;
    DELTA_temp(minInd) = Inf;
    gamma_temp(minInd) = Inf;
    xrange(minInd) = Inf;
    
    clear rect rhomin deltamin D minInd ic
end
clear i cmap RUO_temp DELTA_temp gamma_temp xrange

set(H,'PaperUnits','points','PaperPosition',[200 100 scrsz(4)*1.4 scrsz(4)/2],'PaperPositionMode','manual');
% print(H,'-djpeg',[out,filesep,'Decision_graph.jpg']); clear H

%% assignment
assig = zeros(1,length(RUO));
assig(LabelClust) = LabelClust;

[~,sortInd] = sort(RUO,'descend');
for i=1:length(sortInd)
    if assig(sortInd(i)) == 0
        maxdis = max(dis_matrix(:));
        for j = 1:i-1
            if dis_matrix(sortInd(i),sortInd(j)) < maxdis
                maxdis = dis_matrix(sortInd(i),sortInd(j));
                neigh = sortInd(j);
            end
        end
        assig(sortInd(i)) = assig(neigh);
        clear maxdis j neigh
    end
end
clear i sortInd

% % %% 2D Nonclassical multidimensional scaling
% % H = figure;
% % disp('Performing 2D nonclassical multidimensional scaling');
% % Y1 = mdscale(dis_matrix, 2, 'criterion','metricstress');
% % plot(Y1(:,1),Y1(:,2),'o','MarkerSize',4,'MarkerFaceColor','k','MarkerEdgeColor','k');
% % title ('2D Nonclassical multidimensional scaling','FontSize',10.0);
% % xlabel ('X');
% % ylabel ('Y');
% % 
% % A = zeros(size(Y1,1),2);
% % for i=1:Nclust
% %     nn=0;
% %     for j=1:N
% %         if assig(j) == LabelClust(i)
% %             nn=nn+1;
% %             A(nn,1)=Y1(j,1);
% %             A(nn,2)=Y1(j,2);
% %         end
% %     end
% %     hold on;
% %     plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',4,'MarkerFaceColor',cmap_use(i,:),'MarkerEdgeColor',cmap_use(i,:));
% %     clear nn j
% % end
% % clear i A
% % 
% % % saveas(H,[out,filesep,'Nonclassical_Multidimensional_Scaling.jpg']); clear H
end