邻域网格聚类算法matlab代码及测试数据

%%%%%%%%%%%%%%% Neighborhood Grid Cluster NGC %%%%%%%%%%%%%%%%%%%%%%% % The main idea of the approach is that raw data are divided into grids and % the local density of each grid is simultaneously computed, then starting % with the grid of maximum local density, elements that could be inner ones % of the cluster are gradually added layer-by-layer based on neighborhood relation-ships. % Edit by suozi 2016.05.18 HIT % buaasuozi@126.com % QQ: 379786867 % % The information of the paper and this clustering method is as follows: % @article{Suo2018Neighborhood, % title={Neighborhood grid clustering and its application in fault diagnosis of satellite power system}, % author={Suo, Mingliang and Zhu, Baolong and Zhou, Ding and An, Ruoming and Li, Shunli}, % journal={Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering}, % pages={}, % year={}, % note = {in press, DOI: 10.1177/0954410017751991}, % } % If you use this code, please cite the above paper in your work. Copyright@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [clusterLabel,kindofcluster] = NGC(data,intervalN,delta,bordrate,checknoise) % input: data the data for cluster % intervalN the number of grid, usually it is 80~200 % delta the neighborhood radius, usually it is 0.02~0.08 % bordrate the radius of border, usually is 0.2~0.3 and you could not change it if you don't want. % checknoise check if there are some noise. 0 is not, 1 is yes. 1 is default % output: clusterLabel class label % kindofcluster the number of clusters if nargin < 5 checknoise = 1; end if nargin < 4 checknoise = 1; bordrate = 0.2; end if nargin < 3 checknoise = 1; bordrate = 0.2; delta = 0.05; end if nargin < 2 checknoise = 1; bordrate = 0.2; delta = 0.05; intervalN = 100; end if nargin < 1 disp('At least input one parameter') return end data = myscale(data); % 1��the single point, it is in the boundary of the grid minrho1 = exp(-(sqrt((0.5/intervalN)^2 * size(data,2))/(1/intervalN))^2); % 2��the single point, it is in the center of the grid minrho2 = exp(-0); threshrho = minrho2; % here we choose this one, or you could choose the other one threshgrid = delta*intervalN; % the stop threshold of grid number, which may not be the integer checkdist = threshgrid^2*size(data,2); % the real neighborhood radius, it's 'sqrt' clear minrho1 minrho2 %% the grid partition regionInterval = ones(1,size(data,2)).*(1/intervalN); % the interval of grid tilesize = size(data); tilesize(2) = 1; datasign = floor(data ./ repmat(regionInterval, tilesize)) + 1; % the id of each point clear tilesize %% compute the local density based on Gaussian kernel if size(data,2) == 1 % to do the 1-D situation gaussrho = zeros(ones(1,size(data,2)) .* (intervalN+1), 1); else gaussrho = zeros(ones(1,size(data,2)) .* (intervalN+1)); end idxcell = num2cell(datasign, 1); rhoidx = sub2ind(size(gaussrho), idxcell{:}); clear idxcell for i=1:size(data,1) centerdist = sqrt(sum((data(i,:) - (datasign(i,:) - 0.5).*regionInterval).^2)); % here we use the Euclidean distance gaussrho(rhoidx(i)) = gaussrho(rhoidx(i)) + exp(-(centerdist./regionInterval(1)).^2); end % figure % surf(gaussrho) % grid on % grid minor % xlabel('Dimesion 1');ylabel('Dimension 2');zlabel('\rho') % colorbar % for the valid gaussrho [rhovalue, rho_order] = sort(gaussrho(:),'descend'); validid = find(rhovalue,1,'last'); rho_order_id = rho_order(1:validid); % the index change to subscript rhoorder_sub = zeros(validid,size(data,2)); for i=1:validid % the myind2sub is the function edited by me, the original doesn't suit for my demand. rhoorder_sub(i,:) = myind2sub(size(gaussrho),rho_order_id(i)); end rhoorder_sub_id = 1:1:size(rhoorder_sub,1); %% to expend the neighborhood crntrhoid_v = rhoorder_sub_id; c = 1; % the id of class cluster = zeros(size(gaussrho)); % the sign matrix of cluster while ~isempty(crntrhoid_v) currentCnt = rhoorder_sub(crntrhoid_v(1),:); % from the maximum of the local densities cc = num2cell(currentCnt); % check the noise if checknoise == 1 && gaussrho(cc{:}) <= threshrho noise = rhoorder_sub(crntrhoid_v,:); for i = 1:size(noise,1) cc = num2cell(noise(i,:)); cluster(cc{:}) = -1; % -1 is the noise end clear cc break end cluster(cc{:}) = c; clear cc crntrhoid_v = crntrhoid_v(2:end); bordC = currentCnt; % initialize the border center, the first one is the center of cluster flag = 1; % the flag for going on searching in the current cluster while flag == 1 % this loop is for the current cluster, which will break until there is no element added dist2c = []; % initialize the distance between the center and the elements in the cluster inc_id = []; % it corresponds to dist2c for i = 1:size(bordC,1) crntrhoid_sparse = sparse(crntrhoid_v); [~,~,crntrhoid_v] = find(crntrhoid_sparse); for j = 1:length(crntrhoid_v) dist = sum((rhoorder_sub(crntrhoid_v(j),:) - bordC(i,:)).^2); if dist <= checkdist isin = num2cell(rhoorder_sub(crntrhoid_v(j),:)); cluster(isin{:}) = c; dist2c = [dist2c, sum((rhoorder_sub(crntrhoid_v(j),:) - currentCnt).^2)]; inc_id = [inc_id, crntrhoid_v(j)]; crntrhoid_v(j) = 0; % 0 is delete this one end end end % updating the bordc currentrhoID dist2c if ~isempty(dist2c) % some experiments show that we use max(dist2c) here is much better than other bordid = inc_id(dist2c > bordrate*max(dist2c)); bordC = rhoorder_sub(bordid,:); else % flag = 0; break end % figure % pcolor(cluster) % str = strcat('dist2cN =',num2str(length(dist2c)),32,'bordC =',num2str(size(bordC,1))); % title(str) end c = c + 1; end clear bordC bordid centerdist crntrhoid crntrhoid_v currentCnt dist dist2c... flag i inc_id j validid crntrhoid_sparse isin noise kindofcluster = c-1; str = strcat('totalCluster: ',32,num2str(kindofcluster)); disp(str) str = strcat('totalCluster: ',32,num2str(c-1)); disp(str) %% cluster the original data clusterLabel = zeros(size(data,1),1); for i = 1:size(data,1) crntdata = num2cell(datasign(i,:)); clusterLabel(i) = cluster(crntdata{:}); end clear crntdata i end