一百行写小世界网络和无标度网络

　　曾经觉得能写很长的代码就是厉害的表现，随着学习的深入，对于好的代码有了更加深刻的认识，一个好的代码应该：

　　1、可读性。好的代码不仅能让机器读懂，更要让看代码的人读懂----合理布局，逻辑清晰。

　　2、充分发挥语言的优势，比如接下来的matlab矩阵化编程。

　　3、占用尽可能小的内存，运行尽量快，输出结果尽量容易处理。并在此之间找到平衡。

　　。。。。。

　　需要说明的是，如果没有对复杂网络有基本的了解，下面的东西没有看下去的必要。关于小世界网络和无标度网络，可以参考Watts DJ和Strogatz SH的Collective dynamics of 'small-world' networks(NATURE)以及Barabasi AL和Albert R的 Emergence of scaling in random networks(SCIENCE)。

　　之前，国庆期间曾经用了四百行代码（C语言）实现了无标度网络，之后又用一到两百行代码写完了小世界网络。现在看这些代码，确实受制于C语言语法的限制，所以写出来十分冗杂，而且不易理解。由于最近在学习matlab编程，所以尝试了一下，两个代码加起来竟然不到一百行。而且没有使用matlabBGL之类的工具。而且算法上由于语言变的更加自由，所以算法上也优化了不少（由后面的度分度图就可以看出来）

　　此外，虽然matlab基于JVM,但是只要算法设计的好，还是可以完爆之前的C语言，关于无标度网络的博客中提到计算1000个节点需要35min,经过测试，这个matlab程序处理10000个节点的时间是38min,1000个节点是21.9s！新的程序时间复杂度O(N^2),N是节点数量（由于数据结构没有时间系统学习，时间复杂度可能说错了！）。

　　最后，matlab可以对矩阵进行稀疏化处理，而无标度网络正是一个极其稀疏的网络，所以，借助稀疏化操作，使得matlab也能操作较大的矩阵。

先贴上实现小世界网络的matlab代码：

NodesNum = 1000;
K = 10;
p = 1;
A = sparse(NodesNum, NodesNum);
for i=1:K/2
    A = A + diag(ones(1, length(diag(A, i))), i);
end
A = sparse(A);
for i=1:K/2
    A(i, (NodesNum-K/2+i):NodesNum) = 1;
end
A = A + A';
for i=1:K/2
    P = rand(1,NodesNum);
    P = find(P <= p);
    for j=1:length(P)
        while true
            x = fix(rand()*NodesNum)+1;
            if A(x, P(j)) == 0
                A(x, P(j)) = 1;
                A(P(j), x) = 1;
                break;
            end
        end
        
        if P(j) <= NodesNum-i
            A(NodesNum-i, P(j)) = 0;
            A(P(j), NodesNum-i) = 0;
        else
            A(P(j)-NodesNum+i, P(j)) = 0;
            A(P(j), P(j)-NodesNum+i) = 0;
        end
 
    end
end
clear i j P x
B3 =A^3;
B33 = B3(1:NodesNum+1:end);
B2 = A^2;
B22 = B2(1:NodesNum+1:end);
c = B33./(B22.*(B22-1));
c(isnan(c)) = 0;
C = mean(c);
clear B33 B3 B2 B22

Paths = graphallshortestpaths(tril(A), 'directed', false);
Paths = tril(Paths);
L = sum(sum(Paths)) / nchoosek(NodesNum, 2);

再贴上无标度网络实现的matlab代码

tic
m0 = 10;
m = 5;
NodesNum = 1000;
A = sparse(NodesNum, NodesNum);
A(1:m0,1:m0) = round(rand(m0));
A = tril(A); 
A = A+A';
A = A - diag(diag(A));
for i=m0+1:NodesNum
    Degree = sum(A(1:i-1,1:i-1));
    for j=2:i-1
        Degree(j) = Degree(j) + Degree(j-1);
    end
    LinksNum = 0;
    while LinksNum<m
        link = fix(rand()*Degree(i-1)+1);
        for j=1:i-2
            if link<=Degree(1) && A(i,1)==0
                A(i,1) = 1;     
                A(1,i) = 1;
                LinksNum = LinksNum+1;
            elseif link>Degree(j) && link<=Degree(j+1) && A(i,j+1)==0
                A(i,j+1) = 1;       A(j+1,i) = 1;
                LinksNum = LinksNum+1;
            end
        end
    end
end
Degree = sum(A);
list = unique(Degree);
num = zeros(1,length(list));
for i=1:length(list)
    num(i) = length(find(list(i)==Degree));
end
toc
loglog(list,num ./ sum(num),'.','markersize',20)
xlabel('k'),ylabel('P(k)')

再贴一张小世界网络聚类系数和最短平均路径变化图（最短路径变化有点波动，看来算法还需优化）

最后贴一张10000节点生成的matlab度分布对数坐标图，确实比之前的明显多了（算法改进立竿见影）：

2015-12-14