I need some help in creating a Mathematica program for building and visualizing Watts' "Small Worlds" α-model. This is a variant of the more known (and much more simple) β-model (also known as Watts–Strogatz model), that I think is also implemented in Mathematica, i.e. with the command
RandomGraph[WattsStrogatzGraphDistribution[50, 0.05]]
but I think that the α-model is much more interesting for its intrinsic dynamic characteristics and I'd like to understand it better. The rationale of the α-model is explained in the book:
Duncan J. Watts - Small Worlds: The Dynamics of Networks between Order and Randomness - [1999]
Anyway there's a free preview snippet by the editor at Google Books site that explains it (from page 44 to page 48). It's here (click on page 46) so I don't need to re-explain it here (I just hope the link works also for not "it" (italian) google users).
Till now I've managed to implement the algorithm in Mathematica with the following code:
Module[{n, k, w},
(*
The relevant parameters are;
n: number of nodes;
k: average edges per node (can be decimal);
α: caveman (α\[Rule]0)/solaria (α\[Rule]\
∞) propensity;p: baseline probability for creating edges \
between "distant" nodes (the ones with no mutual friends);
secondary derived parameters are;
M: total number of edges, where M=Round[(k n)/2,0];
iterations: total number of edges to be built after the substrate \
ring initial condition, where iterations=M-n=Round[(k n)/2,0]-n;
*)
Clear@rmat;
Clear@prob;
Clear@cumprob;
Clear@ransamloop;
Clear@matdyn;
Clear@rr;
Clear@αvalue;
Clear@pvalue;
Clear@newnode;
debug = False;
(*debug=True;*)
n = 10;
k = 9;(* it's better if (n k) is even*)
αvalue = 0.1;
pvalue = 0.3;
M = Round[(k n)/2];
iterations = M - n;
activenodes = Range[n];
fuconodes = {};
ransamloop = {};
samfwdcleaned = {};
prob = Table[0, {n}, {n}];
Table[ransamloop =
Flatten@Append[ransamloop, RandomSample[Range[n]]], {s, 1,
Floor[iterations/n] + 1}];
ransamloop = Take[ransamloop, iterations];
matdyn[0] = AdjacencyMatrix[CycleGraph[n]];
mu[i_, j_, matr_] := Module[{count, k},
If[
Or[matr[[i, j]] == 1, i == j], 0,
(*else*)
count = 0;
For[t = 1, t <= n, t++,
If[(matr[[t, i]] == 1) && (matr[[t, j]] == 1), count = count + 1]
]; count
]
];
r[i_, j_, matr_] :=
If[Or[matr[[i, j]] == 1, i == j], 0,(*else*)
If[mu[i, j, matr] >= k,
1,(*else*)
If[mu[i, j, matr] == 0, pvalue,(*else*)
(mu[i, j, matr]/k)^αvalue (1 - pvalue) + pvalue]]];
(* here starts the main loop *)
For[w = 1, w <= iterations, w++,
rmat = Table[Table[r[i, j, matdyn[w - 1]], {j, 1, n}], {i, 1, n}];
fuconodes = {};(*there can be 0, 1,
or 2 nodes completed in the previous step *)
Table[Table[
If[Total[rmat[[i, All]]] == 0,(*
that means that node i is fully connected*)
prob[[i, j]] = 0.;
fuconodes =
DeleteDuplicates[Flatten[Append[fuconodes, i], 1]],(*else*)
prob[[i, j]] = rmat[[i, j]]/Total[rmat[[i, All]]]],
{j, 1, n}],
{i, 1, n}];
(*above if is to avoid that a fully connected node show a prob=0/
0 for all other nodes*)
fuconodes = Flatten[DeleteDuplicates[fuconodes]];
activenodes = DeleteCases[activenodes, Alternatives @@ fuconodes];
todel = Count[ransamloop[[w ;; Length[ransamloop]]],
Alternatives @@ fuconodes];
samfwdcleaned =
DeleteCases[ransamloop[[w ;; Length[ransamloop]]],
Alternatives @@ fuconodes];
toadd = RandomSample[activenodes, todel];
ransamloop = Join[ransamloop[[1 ;; w - 1]], samfwdcleaned, toadd];
rr = RandomReal[];
cumprob =
Table[Table[Sum[prob[[i, j]], {j, 1, u}], {u, 1, n}], {i, 1, n}];
newnode =
FirstPosition[(cumprob[[ransamloop[[w]],
All]](*/.{p -> pvalue, α -> αvalue}*)), _?(# > rr &)][[1]];
If[debug,
Print[{w, "ransamloop->", ransamloop,
"length of ransamloop" -> Length[ransamloop]}]];
If[debug,
Print[Column[{w, (*"rmat->",MatrixForm@rmat,*)"rmat->",
MatrixForm[Round[rmat, 0.01],
TableHeadings -> {Automatic, Automatic}], (*"matdyn[w]->",
MatrixForm@matdyn[w-1],*)"fuconodes->", fuconodes,
"activenodes->", activenodes, "todel->", todel, "toadd->",
toadd, "samfwdcleaned->", samfwdcleaned(*,"rr->",rr*),
"Curr node, rr, newnode->", {ransamloop[[w]], rr, newnode}}(*,
Center*)]]];
(*Print["Curr node, rr, newnode->"{ransamloop[[w]],rr,
newnode}];*)(*matdyn[w]=matdyn[w-1];*) (*not necessary? *)
matdyn[w] =
ReplacePart[matdyn[w - 1],
{{newnode, ransamloop[[w]]}, {ransamloop[[w]], newnode}} -> 1]];
Print["final graph adjacency matrix->",
MatrixForm@matdyn[iterations]];
Print[{"Total number of edges after each step of the main loop->",
1/2 Table[Total[matdyn[f], 2], {f, 1, iterations}]}]]
Then, with the following Manipulate command
Manipulate[
GraphPlot[matdyn[t], VertexLabeling -> True,
VertexCoordinateRules ->
Table[{s -> {Cos[π/2 - (s - 1) (2 π)/Length[matdyn[w]]],
Sin[π/2 - (s - 1) (2 π)/Length[matdyn[w]]]}}, {s, 1,
11}][[All, 1]]], {{t, 0}, 0, iterations, 1,
Appearance -> "Labeled"}]
I can see the network evolution with increasing t (the default parameters n = 10 and k=9 will lead to a fully connected graph for t = 35 given that the starting connected cycle-graph (at t = 0) already has 10 edges and that a fully connected graph with n = 10 has (10*9)/2 = 45 edges).
My problems are:
I've some basic knowledge of Mathematica but I'm absolutely no expert or advanced user.
The program is not optimized at all. It has to do much more calculations than I think are necessary. The processing time is acceptable for n=10 but becomes much longer for higher n (n > 20) and I'd like to check the result for n = 1000 (!) like in Watts' book reported simulations.
I'd like to pack the whole program in a single Manipulate command in which it should be possible to dynamically adjust the relevant parameters (t, k, αvalue, pvalue and even n) and see the corresponding result (my first naive trials to do that led me to infinite loops, errors and weird loop updates of the graph).
Any help?
ParallelTableinstead ofTable? Also consider replaceAppendwith combination ofsowandreap. Also, usually in Mathematica,Doloop is faster thanForloop. $\endgroup$