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

⼀百⾏写⼩世界⽹络和⽆标度⽹络
曾经觉得能写很长的代码就是厉害的表现，随着学习的深⼊，对于好的代码有了更加深刻的认识，⼀个好的代码应该：
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);
5for i=1:K/2
6 A = A + diag(ones(1, length(diag(A, i))), i);
7 end
8 A = sparse(A);
9for i=1:K/2
10 A(i, (NodesNum-K/2+i):NodesNum) = 1;
11 end
12 A = A + A';
13for i=1:K/2
14 P = rand(1,NodesNum);
15 P = find(P <= p);
16for j=1:length(P)
17while true
18 x = fix(rand()*NodesNum)+1;
19if A(x, P(j)) == 0
20 A(x, P(j)) = 1;
21 A(P(j), x) = 1;
22break;
23 end
24 end
25
26if P(j) <= NodesNum-i
27 A(NodesNum-i, P(j)) = 0;
28 A(P(j), NodesNum-i) = 0;
29else
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));
10for i=m0+1:NodesNum
11 Degree = sum(A(1:i-1,1:i-1));
12for j=2:i-1
13 Degree(j) = Degree(j) + Degree(j-1);
14 end
15 LinksNum = 0;
16while LinksNum<m
17 link = fix(rand()*Degree(i-1)+1);
18for j=1:i-2
19if 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));
33for 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。