一百行写小世界网络和无标度网络
曾经觉得能写很长的代码就是厉害的表现，随着学习的深入，对于好的代码有了更加深刻的认识，一个好的代码应该：
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代码：
1 NodesNum = 1000; 2 K = 10; 3 p = 1; 4 A = sparse(NodesNum, NodesNum); 5 for i=1:K/2 6 A = A + diag(ones(1, length(diag(A, i))), i); 7 end 8 A = sparse(A); 9 for i=1:K/2 10 A(i, (NodesNum-K/2+i):NodesNum) = 1; 11 end 12 A = A + A'; 13 for i=1:K/2 14 P = rand(1,NodesNum); 15 P = find(P <= p); 16 for j=1:length(P) 17 while true 18 x = fix(rand()*NodesNum)+1; 19 if A(x, P(j)) == 0 20 A(x, P(j)) = 1; 21 A(P(j), x) = 1; 22 break; 23 end 24 end 25 26 if P(j) <= NodesNum-i 27 A(NodesNum-i, P(j)) = 0; 28 A(P(j), NodesNum-i) = 0; 29 else 30 A(P(j)-NodesNum+i, P(j)) = 0; 31 A(P(j), P(j)-NodesNum+i) = 0; 32 end 33 34 end 35 end 36 clear i j P x 37 B3 =A^3; 38 B33 = B3(1:NodesNum+1:end); 39 B2 = A^2; 40 B22 = B2(1:NodesNum+1:end); 41 c = B33./(B22.*(B22-1)); 42 c(isnan(c)) = 0; 43 C = mean(c); 44 clear B33 B3 B2 B22 45 46 Paths = graphallshortestpaths(tril(A), 'directed', false); 47 Paths = tril(Paths); 48 L = sum(sum(Paths)) / nchoosek(NodesNum, 2);
再贴上无标度网络实现的matlab代码
1 tic 2 m0 = 10; 3 m = 5; 4 NodesNum = 1000; 5 A = sparse(NodesNum, NodesNum); 6 A(1:m0,1:m0) = round(rand(m0)); 7 A = tril(A); 8 A = A+A'; 9 A = A - diag(diag(A)); 10 for i=m0+1:NodesNum 11 Degree = sum(A(1:i-1,1:i-1)); 12 for j=2:i-1 13 Degree(j) = Degree(j) + Degree(j-1); 14 end 15 LinksNum = 0; 16 while LinksNum<m 17 link = fix(rand()*Degree(i-1)+1); 18 for j=1:i-2 19 if link<=Degree(1) && A(i,1)==0 20 A(i,1) = 1; 21 A(1,i) = 1; 22 LinksNum = LinksNum+1; 23 elseif link>Degree(j) && link<=Degree(j+1) && A(i,j+1)==0 24 A(i,j+1) = 1; A(j+1,i) = 1; 25 LinksNum = LinksNum+1; 26 end 27 end 28 end 29 end 30 Degree = sum(A); 31 list = unique(Degree); 32 num = zeros(1,length(list)); 33 for i=1:length(list) 34 num(i) = length(find(list(i)==Degree)); 35 end 36 toc 37 loglog(list,num ./ sum(num),'.','markersize',20) 38 xlabel('k'),ylabel('P(k)')
再贴一张小世界网络聚类系数和最短平均路径变化图（最短路径变化有点波动，看来算法还需优化）
最后贴一张10000节点生成的matlab度分布对数坐标图，确实比之前的明显多了（算法改进立竿见影）：
2015-12-14
