WS小世界模型

clc;close all;clear all;
N=input('请输入网络节点数');
K=input('请输入每个节点的邻居数字(偶数)');
p=input('请输入随机化重连的概率p:');
angle=0:2*pi/N:2*pi-2*pi/N;
%%%%%%%%%%%%%%圆形分布节点%%%%%%%%%%%%%%
for i=1:N
    x(i)=cos(angle(i));
    y(i)=sin(angle(i));
end
figure;
plot(x,y,'bo');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%建立矩阵路线图，不可连接为-1，连接为1，默认为0
route_matrix=zeros(N,N);
for i=1:N
    for j=1:i
         route_matrix(i,j)=-1;%%%我竟不能连我自己
    end
end
%%%%%%%%%%%%%%耦合连接%%%%%%%%%%%%%%
for i=1:N
    for j=i+1:N
        if j-i<=K/2
            route_matrix(i,j)=1;
        elseif i+N-j<=K/2
            route_matrix(i,j)=1;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LineObjects = findall(figure(1),'type','line');
% for i=1:K*N
%     delete(LineObjects(i,1));
% end

%%%%%%%%%%%%%%断边重连%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N-1%%遍历每一个节点，由于是上三角矩阵，所以遍历倒数第二行即可
    b=find(route_matrix(i,:)==1);%%形成每一行有边的数组
    prob=rand(1,1);%%生成一个随机数，判断是否断边
    if length(b)>0%%如果存在边的话
        for j=1:length(b)%%遍历每个边
            if prob>1-p%%需不需要断边
                random_break=rand(1,1);%%生成一个随机数，判断断行点还是断列点
                if random_break>0.5%%断列点
                    a=find(route_matrix(i,:)==0);%%a储存该行可以连的边
                    c=find(route_matrix(:,i)==0);%%c储存该列可以连的边
                    if length(a)+length(c)>0%%如果有空缺可以连边
                        rand_num=randi([1,length(a)+length(c)],1,1);%%生成一个1到空缺边总数的随机整数
                        if rand_num<=length(a)%%在行上有边
                            route_matrix(i,a(rand_num))=1;%%重连
                            route_matrix(i,b(j))=0;%%断边
                        else%%在列上有边
                            route_matrix(c(rand_num-length(a)),i)=1;%%重连
                            route_matrix(i,b(j))=0;%%断边
                        end
                    end
                else%%断行点
                    a=find(route_matrix(b(j),:)==0);
                    c=find(route_matrix(:,b(j))==0);
                    if length(a)+length(c)>0
                        rand_num=randi([1,length(a)+length(c)],1,1);
                        if rand_num<=length(a)
                            route_matrix(b(j),a(rand_num))=1;
                            route_matrix(i,b(j))=0;
                        else
                            route_matrix(c(rand_num-length(a)),b(j))=1;
                            route_matrix(i,b(j))=0;
                        end
                    end
                end
            end
        end
    end
end

                    
%%%%%%%%%%%%%%绘图%%%%%%%%%%%%%%
for i=1:N
    for j=i+1:N
        if route_matrix(i,j)==1
           plot([x(i),x(j)],[y(i),y(j)]);
        end
    end
end
title('WS小世界模型');
hold off;