根据管网信息实现DMA管网自动分区代码案例
利用密度峰值聚类算法代替谱聚类管网分区算法最后一步使用kmeans聚类算法的代码实现
环境要求：（如能直接运行请忽略）
（1）OpenWaterAnalytics/EPANET-Matlab-Toolkit
（2）MATLAB Support for MinGW-w64 C/C++ Compiler
部分运行测试结果图演示如下：
以下代码是在MATLAB2020中运行的，其他版本可能存在某些功能无法使用，可根据要求修改对应的部分。
主文件：
clc
clear
close all
%不同管网相似度矩阵的定义会有不同的管网分区结果，
%% 输入：
sizeCC=30;%点可视化大小
switchColor = 1;%关闭颜色开关
net_name = 'balerma.inp';%管网文件接口
% net_name = 'KL.inp';%
% net_name = 'ky4.inp'
%% 初始定个分区数量，根据管网决策图中密度峰值点与其他点的区别去确定当前分区数量是否合适
PartitionNum =4;%自定义初始分区数量
dc =2;%在多篇DPC聚类算法论文中该参数一般选择为2
c = PartitionNum;%特征向量个数和分区数量有关，论文《基于三维点云的特征分割及其位置关系识别》第三章可推出
is_weightpipe = 1;%相似度矩阵是否增加管道流量权重：1/0
%% 初始信息
d = epanet(net_name);
disp('***************************************************')
disp(['当前管网：',net_name,'分区数量：',num2str(PartitionNum)])
disp('***************************************************')
getConnectivityMatrix = d.getConnectivityMatrix;
getNodeIndex = d.getNodeIndex;
getLinkCount = d.getLinkCount;
getNodeTypeIndex = d.getNodeTypeIndex;
getLinkNameID = d.getLinkNameID;
NodesConnectingLinksIndex = d.getNodesConnectingLinksIndex;
getNodeCoordinates = [d.getNodeCoordinates{1}',d.getNodeCoordinates{2}'];
color_series = {'r','g','b','c','m','y','k',...
[0 0.4470 0.7410],[0.8500 0.3250 0.0980],[0.9290 0.6940 0.1250],[0.4940 0.1840 0.5560],[0.4660 0.6740 0.1880],...
[0.3010 0.7450 0.9330],[0.6350 0.0780 0.1840],[0.4470 0.9410 0.1],[0.3250 0.0980 0.8500],[0.6940 0.1250 0.9290],...
[0.1840 0.5560 0.4940],[0.6740 0.1880 0.4660],[0.7450 0.9330 0.3010],[0.0780 0.1840 0.6350],[0.5 0.4 0.3],...
[0.7 0.2 0.4],[0.8 0.2 0.5],[0.9 0.3 0.1],[0.2 0.1 0.3],[0.9 0.7 0.3],[0.9 0.2 0.3],[0.8 0.7 0.3],[0.1 0.1 0.3],...
[0.8 0.8 0.2],[0.3 0.7 0.7],[0.4 0.7 0.5],[0.3 0.5 0.3],[0.7 0.1 0.9],[0.5 0.7 0.3],[0.4 0.4 0.3],[0.3 0.1 0.3]};%红 绿 蓝 青蓝 品红 黄 黑
coord1 = d.getNodeCoordinates{1};
coord2 = d.getNodeCoordinates{2};
coord = [coord1' coord2'];%物理坐标
%% 水力运算
d.openHydraulicAnalysis;
d.initializeHydraulicAnalysis;
tstep=1;P=[];T_H=[];D=[];H=[];F=[];S=[];HL=[];V=[];
while (tstep > 0) %整点时刻的水力计算
t= d.runHydraulicAnalysis;
if mod(t,3600) ~=0
tstep=d.nextHydraulicAnalysisStep;
continue
end
HL = [HL;d.getLinkHeadloss];
P = [P; d.getNodePressure];
D = [D; d.getNodeActualDemand];
H = [H; d.getNodeHydaulicHead];
S = [S; d.getLinkStatus];
F = [F; d.getLinkFlows];
V = [V;d.getLinkVelocity];%流速
T_H = [T_H; t];
tstep=d.nextHydraulicAnalysisStep;
end
d.closeHydraulicAnalysis;%计算完成
PP = P(1,:)';PPb = PP;
DD = D(1,:)';DDb = DD;
HH = H(1,:)';
FF = F(1,:)';
VV = V(1,:)';
%不同变量往往量纲不同，归一化可以消除量纲对最终结果的影响，是不同变量具有可比性。
FF = FF/sqrt(sum((FF.^2)));
PP = PP/sqrt(sum((PP.^2)));
DD = DD/sqrt(sum((DD.^2)));
HH = HH/sqrt(sum((HH.^2)));
node_features = [PP,DD,HH];
flowMatrix = getConnectivityMatrix;%用于后续根据最短路径确定阀门水表位置
flowMatrix2 = getConnectivityMatrix;
for ii = 1:getLinkCount
link_Weight = abs(FF);
flowMatrix(NodesConnectingLinksIndex(ii,1),NodesConnectingLinksIndex(ii,2)) = 1/link_Weight(ii);
flowMatrix(NodesConnectingLinksIndex(ii,2),NodesConnectingLinksIndex(ii,1)) = 1/link_Weight(ii);
flowMatrix2(NodesConnectingLinksIndex(ii,1),NodesConnectingLinksIndex(ii,2)) = link_Weight(ii);
flowMatrix2(NodesConnectingLinksIndex(ii,2),NodesConnectingLinksIndex(ii,1)) = link_Weight(ii);
end
%% 得到处理结果
%数据特征表示矩阵X -> (node_features = [PP,DD,HH];)
X0 = node_features;
X0 = X0*(X0');
getConnectivityMatrix2 = getConnectivityMatrix + eye(size(getConnectivityMatrix,1));%增加自环
if is_weightpipe
getConnectivityMatrix2 = flowMatrix2.*getConnectivityMatrix2;%%%%%%%%%ddd
disp('增加了流量权重')
end
X0 = getConnectivityMatrix2 .* X0;%构成相似度矩阵
D = diag(sum(X0));
L = D - X0;
Lsym = D^(-1/2) * L * D^(-1/2);
[eigenvector0,eigenvalue0] = eig(Lsym);
% [~,ord] = sort(diag(eigenvalue0),'descend');
[~,ord] = sort(diag(eigenvalue0));
eigenvector = eigenvector0(:,ord(1:c));
Y = normalize(eigenvector,2,'norm');%Y是管网特征空间
%% 节点相似性聚类
%使用DPC算法进行聚类
[idx,C,rho,delta,gamma] = dpc(Y,PartitionNum,dc);%C是密度峰值点坐标；idx是每一点的归属区域
% [idx,C] = kmeans(Y,PartitionNum);%使用Kmeans并没有密度峰值这个概念，不能标记相关内容，且Kmeans算法初始点的选取是随机的，可能该算法产生的多次结果不一致。
%1、确定全局的密度峰值点的标签,双循环找到密度峰值索引peakind
s_C = size(C,1);%密度峰值点数量
s_Y = size(Y,1);
peakind=zeros(1,s_C);
for ii = 1:s_Y
for jj = 1:s_C
yy = Y(ii,:);
cc = C(jj,:);
if yy == cc
peakind(jj) = ii;
end
end
end
disp('密度峰值ID号：')
disp(peakind)
%2、在管网分区结果中可视化密度峰值点
nodeplot = [];colorNodeSet = [];nodeplot2 = [];colorNodeSet2 = [];
for ii = 1:PartitionNum
nodeSetii = d.getNodeNameID(find(idx==ii));
long = length(nodeSetii);
ColornodeSetii = {};
for jj = 1:long
ColornodeSetii(jj)=color_series (ii);
end
nodeplot = [nodeplot nodeSetii];
colorNodeSet = [colorNodeSet ColornodeSetii];
end
% [shortpathLinks,shortpathcolorSet] = plotshortestpath(net_name,peakind);%%%调用函数返回最短路径管道
d.plot('highlightnode',nodeplot,'colornode',colorNodeSet);
hold on;
plotpeak(getNodeCoordinates(peakind,:),150)
hold off;
% d.plot('highlightnode',nodeplot,'colornode',colorNodeSet,'highlightlink',shortpathLinks,'colorlink',shortpathcolorSet);
set(gcf,'unit','normalized','position',[0.275,0.067361111111111,0.37890625,0.815277777777778]);%111111111111111111111
%% 根据密度峰值点开始确定管网阀门水表位置，默认所有边界管道中管道状态关闭的是阀门位置，未关闭的管道是水表位置
%修改储存管网文件
%1.获取边界管道集合
num_border_links = 0;%边界管道数量
for ii = 1:getLinkCount
firstnode = NodesConnectingLinksIndex(ii,1);
secondnode = NodesConnectingLinksIndex(ii,2);
if idx(firstnode) ~= idx(secondnode)%即第ii条管道为边界管道
num_border_links = num_border_links + 1;
border_links(num_border_links) = ii;
end
end
disp(['边界管道数量',num2str(length(border_links))])
disp(['分别为：',num2str(border_links)])
%2.关闭所有边界管道
% d.getLinkInitialStatus(border_links)
status_close = zeros(1,length(border_links));
d.setLinkInitialStatus(border_links, status_close);
% d.getLinkInitialStatus(border_links)
disp('已关闭所有边界管道')
%找到当前分区包括哪些边界管道,cell格式变量
abc1 = NodesConnectingLinksIndex(border_links,:);
for ii = 1:PartitionNum %partition_borderInfo表示属于分区ii的边界管道
temp1 = 0;
for jj = 1:length(border_links)
if (idx(abc1(jj,1)) == ii) || (idx(abc1(jj,2)) == ii)
temp1 = temp1 + 1;
partition_borderInfo{ii}(temp1) = border_links(jj);
end
end
end
%打开各分区某一管道:利用密度峰值点与水源点确定入水口
waterlocal = getNodeIndex(getNodeTypeIndex == 1);%水源位置点索引
Sparse_graph = sparse(flowMatrix);
flag = 1;
while flag
cur_dist = inf;cur_path_index = 0;path_link = [];shortpathLinks = [];shortpathcolorSet = [];tag2 = [];
fixtime = 1;%定义模拟失败时的修改次数
tag5 = 0;%标志:如果为0，则代表还需要修改,需要调整
tag6 = 0;
if tag5 == 0
%%%first
for ii = 1:length(peakind)
for jj = 1:length(waterlocal)
[dist,~,~] = graphshortestpath(Sparse_graph,peakind(ii),waterlocal(jj)); %最短路径图论工具箱
if dist < cur_dist
cur_dist = dist;
cur_path_index = jj;
end
end
%一条最短路径path
[~,path,~] = graphshortestpath(Sparse_graph,peakind(ii),waterlocal(cur_path_index));
%path是点表示的路径，需转换为管道索引并输出
path_length = length(path);
pipe_border_index = partition_borderInfo{ii};%得到第ii个分区的边界管道集合index
for kk = 1:length(pipe_border_index)
abc2 = NodesConnectingLinksIndex(pipe_border_index(kk),:);
if ismember(abc2(1),path) && ismember(abc2(2),path)
needopen_link = pipe_border_index(kk); %要打开的管道
d.setLinkInitialStatus(needopen_link, 1);
disp(['打开了第',num2str(needopen_link),'号边界管道'])
tag2 = [tag2,needopen_link];
end
end
end
%%%second
%储存已关闭的临时计算水力结果的文件
Savefilename = ('fixnet.inp');
d.saveInputFile(Savefilename) ; %存储当前管网结构为inp文件
while tag6==0
% 计算临时存储的文件'fixnet.inp'
f = epanet('fixnet.inp');
f.openHydraulicAnalysis;
f.initializeHydraulicAnalysis;
tstepf=1;Pf=[];
while (tstepf > 0) %整点时刻的水力计算
t= f.runHydraulicAnalysis;
if mod(t,3600) ~=0
tstepf=f.nextHydraulicAnalysisStep;
continue
end
Pf = [Pf; f.getNodePressure];
tstepf=f.nextHydraulicAnalysisStep;
end
f.closeHydraulicAnalysis;%计算完成
PPf = Pf(1,:)' ; %分区后
%判断是否模拟成功
tag4 = sum(PPf)-sum(abs(PPf));
if abs(tag4) < 1
disp('***************************************************')
disp(['当前调整模拟成功，共经历',num2str(fixtime),'次修改'])
disp('***************************************************')
flag = 0;
tag6 = 1;
else
disp(['还可开启的管道:',getLinkNameID{setdiff(border_links,tag2)}])%无用
disp('??????????????????????????????????????????')
disp('模拟出现异常，需要增加开启的管道')
disp(['修改次数为：',num2str(fixtime)])
%如果模拟失败，打开一条流量大的边界管道
fmax = 0;
abc3 = setdiff(border_links,tag2);%还可以再打开的管道???????此处要更新
for uu = 1:length(abc3)%判断流量大小
if abs(FF(abc3(uu))) > fmax
fmax = abs(FF(abc3(uu)));
fmaxpipe = abc3(uu);
disp(['最大流量管道：',getLinkNameID(fmaxpipe)])
end
end
%排除abc3中的打开的管道
tag2 = [tag2 , fmaxpipe];
f.setLinkInitialStatus(fmaxpipe,1);
disp(['开启了第',num2str(fmaxpipe),'号',getLinkNameID{fmaxpipe},'管道'])
f.saveInputFile(Savefilename) ; %存储当前管网结构为inp文件
end
fixtime = fixtime + 1;
end
else
disp('----------------')
disp(['不需要调整了，已调整：',num2str(fixtime)])
break;
end
end
disp('------------------------------------------------------------------')
%3 计算压力相关数据(分区前后)
PPb;DDb;%可显示分区前数据
%分区后；
f = epanet('fixnet.inp');%获取已修改完成的可运行的管网的数据
f.openHydraulicAnalysis;
f.initializeHydraulicAnalysis;
tstepf=1;Pf=[];Df=[];
while (tstepf > 0) %整点时刻的水力计算
t= f.runHydraulicAnalysis;
if mod(t,3600) ~=0
tstepf=f.nextHydraulicAnalysisStep;
continue
end
Df = [Df; f.getNodeActualDemand];
Pf = [Pf; f.getNodePressure];
tstepf=f.nextHydraulicAnalysisStep;
end
f.closeHydraulicAnalysis;%计算完成
PPf = Pf(1,:)' ; %分区后
DDf = Df(1,:)' ; %分区后
PPb;%分区前
std0 = zeros(1,PartitionNum);std000 = zeros(1,PartitionNum);
%还得计算各分区最高最低压力以及压力标准差/边界管道数量。
disp('---------------------------分区后显示数据开始---------------------------------')
for ii = 1:PartitionNum
temp2 = PPf(idx == ii);%分区后当前分区的压力
temp4 = DDf(idx == ii);%分区后当前分区的需水量
PartitionNodenum = length(temp2);%%%%%%%%%%%%%%%其他版本未改过
temp2 = temp2(temp2>0.01);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%其他版本未改过
DDSource = temp4(temp4<-0.01);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%其他版本未改过
avg_PPf = mean(temp2);
max_F = max(temp2);%%%%%%%%%%%%%%%其他版本未改过
min_F = min(temp2);%%%%%%%%%%%%%%%其他版本未改过
disp(['| ','DMA',num2str(ii),'分区：','边界管道数量->',num2str(length(partition_borderInfo{ii}))])
disp(['| ','DMA',num2str(ii),'分区：','平均压力->',num2str(avg_PPf),' ','节点数量->',num2str(PartitionNodenum)])%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%其他版本未改过
disp(['| ','DMA',num2str(ii),'分区：','最高压力以及最低压力->',num2str(max_F),' ',num2str(min_F)])%%%%%%%%%%%%%%%其他版本未改过
diff1 = (temp2 - avg_PPf).^2;
std1 = sqrt(sum(diff1)/length(temp2));
disp(['| ','DMA',num2str(ii),'分区：','分区节点压力标准差：',num2str(std1)])%%%%%%%%%%%%%%%其他版本未改过
std0s(ii) = std1;
disp(' +++++++++++++++++++++++ ')
end
disp('---------------------------分区后显示数据结尾---------------------------------')
std0 = mean(std0s);
disp('-----------------------------------------')
disp(['分区后各分区的压力标准差的平均值',num2str(std0)])
disp('-----------------------------------------')
std00 = sqrt(mean((PPb - mean(PPb)).^2));
disp(['分区前整块分区的压力标准差',num2str(std00)])
%% 可视化密度峰值点分布
figure;
if switchColor
for ii = 1:PartitionNum
locidx = find(idx == ii);
hold on ;
% colo = color_series{3};
colo = color_series{ii};%可视化版
scatter(rho(locidx),delta(locidx),'filled','MarkerFaceColor',colo);
hold off;
end
hold on ;
scatter(rho(peakind),delta(peakind),120,'kp','LineWidth',1);%md
hold off;
else
scatter(rho,delta,'filled','b');
end
xlabel('局部密度');
ylabel('相对距离');
title('管网决策图');
figure;%特征空间
YY = Y;
for ii = 1:PartitionNum
locidx = find(idx == ii);
xx = YY(locidx,1);
yy = YY(locidx,2);
zz = YY(locidx,3);
% colo = color_series{ii};%可视化版
colo = color_series{3};
scatter3(xx,yy,zz,sizeCC,'filled','MarkerFaceColor',colo);
hold on
end
hold off
%% 结构模块度
WW0 = abs(X0);%0.91699；.....0.89544----0.97957+++++++++++++0.88131
WW = sum(sum(WW0))/2;
WFQ_list = zeros(1,PartitionNum);
for jjj = 1:PartitionNum
WWC = 0;WSC = 0;
for iii = 1:getLinkCount
if idx(NodesConnectingLinksIndex(iii,1)) == idx(NodesConnectingLinksIndex(iii,2)) && (idx(NodesConnectingLinksIndex(iii,1))==jjj || idx(NodesConnectingLinksIndex(iii,2))==jjj)
WWCpoint = WW0(NodesConnectingLinksIndex(iii,1),NodesConnectingLinksIndex(iii,2));
WWC = WWCpoint + WWC;
end
if idx(NodesConnectingLinksIndex(iii,1))~= idx(NodesConnectingLinksIndex(ii,2)) && (idx(NodesConnectingLinksIndex(iii,1))==jjj || idx(NodesConnectingLinksIndex(iii,2))==jjj)
WSCpoint = WW0(NodesConnectingLinksIndex(iii,1),NodesConnectingLinksIndex(iii,2));
WSC = WSCpoint + WSC;
end
end
WFQ_list(jjj) = WWC/WW-(WSC/(2*WW))^2;
end
F_Q = sum(WFQ_list);
disp(['结构模块度：',num2str(F_Q),'---------------------------'])
%% 聚类有效性指标函数及其改进
% 设Si表示类i中各点的分离度，表示如下：
Si = zeros(1,PartitionNum);
for ii = 1:length(Si)
Si_Xi = getNodeIndex(idx == ii);
for jj = 1:length(Si_Xi)
Si(ii) = (sum(sum(abs(Y(Si_Xi,:)-Y(peakind(ii),:)),2))) / length(Si_Xi); %######################
end
end
%簇间分离度Sepij:两簇间最远边界点的距离，
Sepij = zeros(PartitionNum,PartitionNum);
for ii = 1:PartitionNum
Sepij_Y1 = Y(getNodeIndex(idx == ii),:);
for jj = 1:PartitionNum
Sepij_Y2 = Y(getNodeIndex(idx == jj),:);
dij_dist = pdist([Sepij_Y1;Sepij_Y2],'cityblock');
Sepij2dist = squareform(dij_dist);
% [Sepij_temp1,Sepij_temp2] = max(Sepij);%Sepij_temp1每一列的最大值，Sepij_temp2每一列最大值的位置
% [Sepij_temp3,Sepij_temp4] = max(Sepij_temp1);%Sepij_temp3最大值，Sepij_temp4最大值位置
%最大值坐标（Sepij_temp4，Sepij_temp2）
if ii ~= jj
Sepij(ii,jj) = max(max(Sepij2dist));
end
end
end
%相似度Rij
Rij = zeros(PartitionNum,PartitionNum);
for ii = 1:PartitionNum
for jj = 1:PartitionNum
if ii~=jj
Rij(ii,jj) = (Si(ii)+Si(jj))/Sepij(ii,jj);
end
end
end
IDBI = max(max(Rij))%值越小越好
% IDBI = mean(max(Rij)) %########################
%%
if is_weightpipe
disp('含有流量权重信息的')
end
%节点均衡性
PART=zeros(1,PartitionNum);
for i = 1:PartitionNum
PART(i) = sum(idx == i);
end
DIFF = (PART - mean(PART)).^2;
III2=sqrt(mean(DIFF));
disp(['节点均衡性',num2str(III2)])
disp(['边界管道数量',num2str(length(border_links))])
%% 分区前后结果对比
%画分区前后压力差
sizesss = 2;
fontsizesss = 15;
figure;
plot(PPb(PPb>0),'r','LineWidth',sizesss)
hold on
plot(PPf(PPf>0),'b','LineWidth',sizesss)
hold off
legend('分区前','分区后','FontSize',fontsizesss,'LineWidth',sizesss)
xlabel('管网节点索引号','LineWidth',sizesss);
ylabel('节点压力(m)','LineWidth',sizesss);
set(gca,'FontSize',fontsizesss,'LineWidth',sizesss);%坐标轴的字体及大小设置
figure;
hold on
plot(PPf(PPf>0)-PPb(PPb>0),'LineWidth',sizesss)
% plot(PPb(PPb>0)-PPf(PPf>0))
hold off
xlabel('管网节点索引号');
ylabel('节点压力差(m)');
set(gca,'FontSize',fontsizesss,'LineWidth',sizesss);%坐标轴的字体及大小设置
%% 函数调用
function plotpeak(peakcoord,sizep)
x = peakcoord(:,1);
y = peakcoord(:,2);
scatter(x,y,sizep,'rp','fill')
end
dpc.m文件：
function [idx,C,rho,delta,gamma] = dpc(Y,K,p)
% 从文件中读取数据
data_load=Y; %%%可更换为自己的数据
mdist=pdist(data_load); %两两行之间距离
A=tril(ones(size(data_load,1)))-eye(size(data_load,1)); %%%得到下三角矩阵
[x,y]=find(A~=0); %%%得到可以标识任何点的坐标x和y
% Column 1: id of element i, Column 2: id of element j', Column 3: dist(i,j)'
xx=[x y mdist'];
ND=max(xx(:,2));
NL=max(xx(:,1));
if (NL>ND)
ND=NL; %% 确保 ND 取为第一二列最大值中的较大者，并将其作为数据点总数
end
N=size(xx,1); %% xx 第一个维度的长度，相当于文件的行数（即距离的总个数）
% 初始化为零
dist=zeros(ND,ND);
% 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
% 1.确定 dc
percent=p; %%%原版是percent=2.0；现在我改为4.0；
fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent);
%%%取所有排序后距离总数量的百分之四的位置作为dc的取值
position=round(N*percent/100); %% round 是一个四舍五入函数
sda=sort(xx(:,3)); %% 对所有距离值作升序排列
dc=sda(position);
% 2.计算局部密度 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.txt')
disp('column 1:Density')
disp('column 2:Delta')
fid = fopen('DECISION_GRAPH.txt', 'w');
for i=1:ND
fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));
end
%%%上面是得到了rho和delta
%%%下面一小段是本人对rho和delta进行归一化，防止数量级不统一带来的影响，先注释
rho = normalize(rho,'range');
delta = normalize(delta,'range');
% % % % % % % 数据归一化第二种方法
% % % % % % rho=mapminmax(rho,0,1);
% % % % % % delta=mapminmax(delta,0,1);
% ind 和 gamma 在后面并没有用到
for i=1:ND
ind(i)=i;
gamma(i)=rho(i)*delta(i); %%%该乘积值可用于‘自动确定聚类数目’，排序后可取前几个
end
% figure();
% hold on ;
% scatter(rho,delta);
% hold off;
% figure();
% hold on ;
% plot(gamma,'bo');
% hold off;
%对gamma进行排序，然后选择突变（跳变）的那些。在这里只选择前 7个最大的
[gamma_sorted,gamma_order]=sort(gamma,'descend');
cccc=gamma_order(1:K); %% KK决定聚类中心个数
C = Y(cccc,:);
% 初始化 cluster 个数
NCLUST=0;
% cl 为归属标志数组，cl(i)=j 表示第 i 号数据点归属于第 j 个 cluster
for i=1:ND
cl(i)=-1;
end
% 直接利用cccc数组获取 聚类中心点，最大的在第一个,其它的非聚类中心点都是-1
[~,NCLUST]=size(cccc);
for i=1:NCLUST
cl(cccc(i))=i;
icl(i)=cccc(i);
end
% 将其他数据点归类 (assignation)
for i=1:ND
if (cl(ordrho(i))==-1) %%%意思是：如果cl（i）= -1，则代表还未归类
cl(ordrho(i))=cl(nneigh(ordrho(i))); %%%nneight代表与i点离得最近的点的编号
end
end
% 由于是按照 rho 值从大到小的顺序遍历,循环结束后, cl 应该都变成正的值了.
idx=cl;
end
1241
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
