-
Notifications
You must be signed in to change notification settings - Fork 77
/
DoddsWattsSabel.m
executable file
·80 lines (64 loc) · 2.47 KB
/
DoddsWattsSabel.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
% Add random cross-links on top of a perfect hierarchy.
% Non-backbone edges are added with probability
% P(i,j)=e^(-Dij/lambda)*e^(-xij/ksi),
% where Dij is the level of the lowest common ancestor
% and xij is the "organizational" distance
% Source: Dodds, Watts, Sabel, "Information exchange and the
% robustness of organizational networks", PNAS 100 (21): 12516-12521
%
% INPUTs: number of nodes (N), tree branch factor (b),
% m - number of additional edges to add,
% parameters lambda (lam) and ksi in [0,inf)
% OUTPUTs: adjacency matrix of randomized hierarchy, NxN
%
% Other routines used: edgeL2adj.m, canonicalNets.m, dijkstra.m
% Last updated: November 23 2012
function adj = DoddsWattsSabel(N, b, m, lam, ksi)
% construct a tree with N nodes and branch factor b
adj0 = edgeL2adj(canonicalNets(N, 'tree', b)); % backbone adjacency matrix
adj = adj0;
edges = 0;
while edges < m
% pick two nodes at random
ind1 = randi(N);
ind2 = randi(N);
% if same node or already connected, keep going
if ind1 == ind2 || adj(ind1, ind2) + adj(ind2, ind1) > 0;
continue;
end
% find di,dj and Dij
[d1, path1] = dijkstra(adj0, ind1, 1); % adjacency, source, target
[d2, path2] = dijkstra(adj0, ind2, 1);
for p = 1:length(path1)
p1 = path1(p);
p2 = find(path2 == p1);
if length(p2) > 0% if p1 in path2
di = p - 1; dj = p2 - 1; % di+dj is the the distance from ind1
% to ind2 along the backbone hierarchy
Dij = length(path1(p:length(path1))) - 1;
% the level of the highest common node on the path to the root
break
end
end
xij = sqrt(di^2 + dj^2 - 2);
% connect ind1 and ind2 with prob. e^(-Dij/lam)e^(-xij/ksi))
if rand < exp(-Dij / lam) * exp(-xij / ksi)
adj(ind1, ind2) = 1;
adj(ind2, ind1) = 1;
edges = edges + 1;
end
end
%!test
%! for x=1:40
%! randint = randi(50)+2;
%! m = randi(round(randint/4));
%! adj = DoddsWattsSabel(randint, 2, m, 10*rand, 10*rand);
%! assert(numEdges(adj), m+randint-1)
%! assert(isTree(adj), false)
%! assert(isConnected(adj), true)
%! assert(isSimple(adj), true)
%! end
%!demo
%! adj = DoddsWattsSabel(50, 3, 10, 15, 15);
%! assert(numEdges(adj), 59)
%! isSimple(adj)