clear all
close all
disp('The only input needed is a distance matrix file')
disp('The format of this file should be: ')
disp('Column 1: id of element i')
disp('Column 2: id of element j')
disp('Column 3: dist(i,j)') %% 從文件中讀取數據
mdist=input('name of the distance matrix file (with single quotes)?\n');
disp('Reading input distance matrix')
xx=load(mdist);
ND=max(xx(:,2));
NL=max(xx(:,1));
if (NL>ND) ND=NL; %% 確保 DN 取為第一二列最大值中的較大者，并將其作為數據點總數
end N=size(xx,1); %% xx 第一個維度的長度，相當于文件的行數（即距離的總個數） %% 初始化為零
for i=1:ND for j=1:ND dist(i,j)=0; end
end %% 利用 xx 為 dist 數組賦值，注意輸入只存了 0.5*DN(DN-1) 個值，這里將其補成了滿矩陣
%% 這里不考慮對角線元素
for i=1:N ii=xx(i,1); jj=xx(i,2); dist(ii,jj)=xx(i,3); dist(jj,ii)=xx(i,3);
end %% 確定 dc percent=2.0;
fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent); position=round(N*percent/100); %% round 是一個四舍五入函數
sda=sort(xx(:,3)); %% 對所有距離值作升序排列
dc=sda(position); %% 計算局部密度 rho (利用 Gaussian 核) fprintf('Computing Rho with gaussian kernel of radius: %12.6f\n', dc); %% 將每個數據點的 rho 值初始化為零
for i=1:ND rho(i)=0.;
end % Gaussian kernel
for i=1:ND-1 for j=i+1:ND rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc)); rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc)); end
end % "Cut off" kernel
%for i=1:ND-1
% for j=i+1:ND
% if (dist(i,j)<dc)
% rho(i)=rho(i)+1.;
% rho(j)=rho(j)+1.;
% end
% end
%end %% 先求矩陣列最大值，再求最大值，最后得到所有距離值中的最大值
maxd=max(max(dist)); %% 將 rho 按降序排列，ordrho 保持序
[rho_sorted,ordrho]=sort(rho,'descend'); %% 處理 rho 值最大的數據點
delta(ordrho(1))=-1.;
nneigh(ordrho(1))=0; %% 生成 delta 和 nneigh 數組
for ii=2:ND delta(ordrho(ii))=maxd; for jj=1:ii-1 if(dist(ordrho(ii),ordrho(jj))<delta(ordrho(ii))) delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj)); nneigh(ordrho(ii))=ordrho(jj); %% 記錄 rho 值更大的數據點中與 ordrho(ii) 距離最近的點的編號 ordrho(jj) end end
end %% 生成 rho 值最大數據點的 delta 值
delta(ordrho(1))=max(delta(:)); %% 決策圖 disp('Generated file:DECISION GRAPH')
disp('column 1:Density')
disp('column 2:Delta') fid = fopen('DECISION_GRAPH', 'w');
for i=1:ND fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));
end %% 選擇一個圍住類中心的矩形
disp('Select a rectangle enclosing cluster centers') %% 每臺計算機，句柄的根對象只有一個，就是屏幕，它的句柄總是 0
%% >> scrsz = get(0,'ScreenSize')
%% scrsz =
%% 1 1 1280 800
%% 1280 和 800 就是你設置的計算機的分辨率，scrsz(4) 就是 800，scrsz(3) 就是 1280
scrsz = get(0,'ScreenSize'); %% 人為指定一個位置，感覺就沒有那么 auto 了 :-)
figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]); %% ind 和 gamma 在后面并沒有用到
for i=1:ND ind(i)=i; gamma(i)=rho(i)*delta(i);
end %% 利用 rho 和 delta 畫出一個所謂的“決策圖” subplot(2,1,1)
tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');
title ('Decision Graph','FontSize',15.0)
xlabel ('\rho')
ylabel ('\delta') subplot(2,1,1)
rect = getrect(1);
%% getrect 從圖中用鼠標截取一個矩形區域， rect 中存放的是
%% 矩形左下角的坐標 (x,y) 以及所截矩形的寬度和高度
rhomin=rect(1);
deltamin=rect(2); %% 作者承認這是個 error，已由 4 改為 2 了! %% 初始化 cluster 個數
NCLUST=0; %% cl 為歸屬標志數組，cl(i)=j 表示第 i 號數據點歸屬于第 j 個 cluster
%% 先統一將 cl 初始化為 -1
for i=1:ND cl(i)=-1;
end %% 在矩形區域內統計數據點（即聚類中心）的個數
for i=1:ND if ( (rho(i)>rhomin) && (delta(i)>deltamin)) NCLUST=NCLUST+1; cl(i)=NCLUST; %% 第 i 號數據點屬于第 NCLUST 個 cluster icl(NCLUST)=i;%% 逆映射,第 NCLUST 個 cluster 的中心為第 i 號數據點 end
end fprintf('NUMBER OF CLUSTERS: %i \n', NCLUST); disp('Performing assignation') %% 將其他數據點歸類 (assignation)
for i=1:ND if (cl(ordrho(i))==-1) cl(ordrho(i))=cl(nneigh(ordrho(i))); end
end
%% 由于是按照 rho 值從大到小的順序遍歷,循環結束后, cl 應該都變成正的值了. %% 處理光暈點，halo這段代碼應該移到 if (NCLUST>1) 內去比較好吧
for i=1:ND halo(i)=cl(i);
end if (NCLUST>1) % 初始化數組 bord_rho 為 0,每個 cluster 定義一個 bord_rho 值 for i=1:NCLUST bord_rho(i)=0.; end % 獲取每一個 cluster 中平均密度的一個界 bord_rho for i=1:ND-1 for j=i+1:ND %% 距離足夠小但不屬于同一個 cluster 的 i 和 j if ((cl(i)~=cl(j))&& (dist(i,j)<=dc)) rho_aver=(rho(i)+rho(j))/2.; %% 取 i,j 兩點的平均局部密度 if (rho_aver>bord_rho(cl(i))) bord_rho(cl(i))=rho_aver; end if (rho_aver>bord_rho(cl(j))) bord_rho(cl(j))=rho_aver; end end end end %% halo 值為 0 表示為 outlier for i=1:ND if (rho(i)<bord_rho(cl(i))) halo(i)=0; end end end %% 逐一處理每個 cluster
for i=1:NCLUST nc=0; %% 用于累計當前 cluster 中數據點的個數 nh=0; %% 用于累計當前 cluster 中核心數據點的個數 for j=1:ND if (cl(j)==i) nc=nc+1; end if (halo(j)==i) nh=nh+1; end end fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(i),nc,nh,nc-nh); end cmap=colormap;
for i=1:NCLUST ic=int8((i*64.)/(NCLUST*1.)); subplot(2,1,1) hold on plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
end
subplot(2,1,2)
disp('Performing 2D nonclassical multidimensional scaling')
Y1 = mdscale(dist, 2, 'criterion','metricstress');
plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','MarkerEdgeColor','k');
title ('2D Nonclassical multidimensional scaling','FontSize',15.0)
xlabel ('X')
ylabel ('Y')
for i=1:ND A(i,1)=0.; A(i,2)=0.;
end
for i=1:NCLUST nn=0; ic=int8((i*64.)/(NCLUST*1.)); for j=1:ND if (halo(j)==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',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
end %for i=1:ND
% if (halo(i)>0)
% ic=int8((halo(i)*64.)/(NCLUST*1.));
% hold on
% plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
% end
%end
faa = fopen('CLUSTER_ASSIGNATION', 'w');
disp('Generated file:CLUSTER_ASSIGNATION')
disp('column 1:element id')
disp('column 2:cluster assignation without halo control')
disp('column 3:cluster assignation with halo control')
for i=1:ND fprintf(faa, '%i %i %i\n',i,cl(i),halo(i));
end
