2
$\begingroup$

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:

  1. I've some basic knowledge of Mathematica but I'm absolutely no expert or advanced user.

  2. 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.

  3. 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?

$\endgroup$
  • 1
    $\begingroup$ I am not an expert by any means but have you tried ParallelTable instead of Table? Also consider replace Append with combination of sow and reap. Also, usually in Mathematica, Do loop is faster than For loop. $\endgroup$ – Mahdi Oct 6 '14 at 5:19
  • $\begingroup$ Thanks @Mahdi for your helpful suggestions. I'll try them. $\endgroup$ – Luca M Oct 6 '14 at 14:25
0
$\begingroup$

I've managed to do some optimization of the code, using matrix operations instead of loops (or tables) whenever possible and changing the main loop from a For to a Do as suggested. The ParallelTable command (also suggested) seems not to work inside a Manipulate.

The computation time has reduced by a factor of 10 (more or less) and it's now possible to build a graph with up to 50 nodes (about 8s wait).

The heavier computations are made in the external Manipulate: changing those parameters triggers the whole graph to be recalculated. The t parameter in the nested Manipulate, instead, shows the steps through which the graph is built till it's final degree and doesn't require a recalculation of the whole graph.

I think that I cannot do more than that (also given my limited knowledge of Mathematica) so I'm posting here below the code I finally came up with. Any suggestion for improvement of the code is welcome.

Manipulate[
 (*
 The relevant parameters are;
 n: number of nodes;
 k: final average edges per node (graph average degree - can be \
decimal);
 \[Alpha]value: caveman (\[Alpha]\[Rule]0)/solaria (\[Alpha]\[Rule]\
\[Infinity]) propensity;pvalue: baseline probability for creating \
edges between "distant" nodes (the ones with no mutual friends);
 Above parameters are set with controls in the external Manipulate. \
Changing one of above parameters triggers the recalculation of the \
graph;
 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;
  *)
 debug = False; (* with debug=True some partial data will be printed \
out during the main cycle to check the process;*)
 (*debug=True;*) (* with debug=True some partial data will be printed \
out during the main cycle to check the process;*)
 Clear@rmat;(*ramt: propensity matrix to establish a connection from \
node i to node j *)
 Clear@mmat;(*"mutual friends" matrix. Subscript[m, ij] is the number \
of mutual friends between node i and node j. 
 If i and j are directly connected than Subscript[m, ij]=0 *)
 Clear@matdyn; (* Indexed adjacency matrix of the graph during the \
mainloop. If matdyn[w][[i,j]]==0 then node i and node j are not \
linked at step w. If matdyn[w][[i,j]]\[Equal]1 then node i and node j \
are linked at step w *)
 Clear@newnode; (* node to be connected to the currently selected \
node (chosen in ransamloop) *)
 If[k > n - 1, k = n - 1];
M = Round[(k n)/2]; (* final number of edges in the graph *)
 iterations = 
  M - n; (* number of the main loop iterations needed to reach M \
edges after the initial graph (ring with k=2 and n edges*)
 activenodes = 
  Range[n]; (* list of the nodes that are not (yet) fully connected - \
updated during the main loop *)
 fuconodes = {}; (*list of the nodes that have become fully completed \
(that is linked to all the other nodes - updated during the main loop \
*)
 onesmat = 
  Table[1, {i, 1, n}, {j, 1, 
    n}]; (* a matrix filled with "1"s in every cell - useful for some \
matrix operation *)
 ransamloop = {};(* each node will be selected in random order, but \
once a node has been allowed to choose a new neighbour (newnode), it \
may not choose again until all other nodes have taken their turn once \
in turn. The choice will be a random choice weighted by the rmat \
propensity values. ransamloop is the sequence of nodes allowed to \
choose their new partners (newnode) and might be updated in each loop \
cycle whenever some node becomes completed after a new link is \
established *)
 samfwdcleaned = {}; (* forward part of ransamloop cleaned from fully \
completed nodes *)

 (*construction of the node sequence to reach M edges, that is M-n \
iterations *)
 Table[ransamloop = 
   Flatten@Append[ransamloop, RandomSample[Range[n]]], {s, 1, 
   Floor[iterations/n] + 1}];
 ransamloop = Take[ransamloop, iterations];
 matdyn[0] = 
  AdjacencyMatrix[
   CycleGraph[n]]; (* the starting graph is a simple ring graph *)
 loopnodecomplete = 
  False; (*boolean to track if the node chosen from ransamloop during \
the cycle has become a completed node *)
 connodecomplete = 
  False; (*boolean to track if the node linked to the node chosen \
from ransamloop during the cycle has become a completed node *)
 regularizer[matr_] := 
  1 - (matr + 
     IdentityMatrix[
      Length[matr]]); (*This function builds a matrix that, if \
multiplied element-wise by an adjacency matrix or by the rmat matrix, \
sets the diagonal elements and the already connected nodes elements \
to 0 to exclude the possibility of "self connected nodes" and \
"already connected nodes" to be re-considered for new connections *)

 (* here starts the main "Do" loop (from w=1 to iterations) *)
 If[debug, 
  Print[Column[{"======= INITIAL CONDITIONS: =======================",
      Row[{"ransamloop->", ransamloop, 
       " - length of ransamloop:" -> Length[ransamloop]}]}]]];
 (* begin of main loop *)
 Do[
  kprog = 2 + 2/n (w - 1);(* average graph degree at step w-1; 
   unsure if to use the final (target) k instead, 
   as Watts' book hints*)
   progr = Round[N[w/iterations 100], 1];
   loopnodecomplete = False;
   connodecomplete = False;
   todel1 = 0; todel2 = 0; samfwdcleaned = {}; toadd1 = {}; 
   toadd2 = {};
   mmat = 
    Table[Table[
       Total[matdyn[w - 1][[j, All]]* matdyn[w - 1][[i, All]]], {j, 1,
         n}], {i, 1, n}]*regularizer[matdyn[w - 1]];
   rmat = 
    If[\[Alpha]value == 0, 
     Table[If[mmat[[i, j]] == 0, pvalue, 1], {i, 1, n}, {j, 1, n}]*
      regularizer[matdyn[w - 1]],(*else: \[Alpha]value<>
     0 *)((MapThread[Min, {mmat/kprog, onesmat}, 
           2])^\[Alpha]value (1 - pvalue) + pvalue)*
      regularizer[matdyn[w - 1]]];
   newnode = RandomChoice[rmat[[ransamloop[[w]], All]] -> Range[n]];
   matdyn[w] = 
    ReplacePart[
     matdyn[w - 
       1], {{newnode, ransamloop[[w]]}, {ransamloop[[w]], newnode}} ->
       1];(* in the following section there is the possible cleanup \
of fully completed nodes from ransamloop and the corresponding \
updates of fuconodes and activenodes *)
   If[Total[matdyn[w][[ransamloop[[w]], All]]] == n - 1,(* 
    that would mean that node "ransamloop[[w]]" is complete *)
    AppendTo[fuconodes, ransamloop[[w]]]; loopnodecomplete = True(*,
    else: do nothing*)]; 
   If[Total[matdyn[w][[newnode, All]]] == n - 1,(* 
    that would mean that node "newnode" is complete *)
    AppendTo[fuconodes, newnode];
    connodecomplete = True(*,else: do nothing*)];.
   If[loopnodecomplete,(*update activenodes and ransamloop *)
    activenodes = DeleteCases[activenodes, ransamloop[[w]]];
    todel1 = 
     Count[ransamloop[[w + 1 ;; Length[ransamloop]]], ransamloop[[w]]];
    samfwdcleaned = 
     DeleteCases[ransamloop[[w + 1 ;; Length[ransamloop]]], 
      ransamloop[[w]]];
    toadd1 = RandomSample[activenodes, todel1];
    ransamloop = 
     Flatten@Join[ransamloop[[1 ;; w]], samfwdcleaned, toadd1](*,
    else: do nothing*)];
  If[connodecomplete,(*update activenodes and ransamloop *)
   activenodes = DeleteCases[activenodes, newnode];
   todel2 = Count[ransamloop[[w + 1 ;; Length[ransamloop]]], newnode];
   samfwdcleaned = 
    DeleteCases[ransamloop[[w + 1 ;; Length[ransamloop]]], newnode];
   toadd2 = RandomSample[activenodes, todel2];
   ransamloop = Join[ransamloop[[1 ;; w]], samfwdcleaned, toadd2](*,
   else: do nothing*)];

  (* loop print section in debug mode BEGIN ========================= \
*)
  If[debug, 
   Print[Column[{Print[
       "==================================================="], 
      Row[{"w = " <> ToString[w], 
        "current ransamloop->" <> ToString@ransamloop[[w]], 
        "length of ransamloop" -> Length[ransamloop]}, "; "], 
      "kprog = " <> ToString[N[kprog]], 
      Row[{"matdyn[w-1]->", 
        MatrixForm[matdyn[w - 1], 
         TableHeadings -> {Automatic, Automatic}], "mmat->", 
        MatrixForm[mmat, TableHeadings -> {Automatic, Automatic}], 
        "  rmat->", 
        MatrixForm[Round[rmat, 0.01], 
         TableHeadings -> {Automatic, Automatic}](*,
        "  probmat->",*)(*MatrixForm[Round[probmat,0.001],
        TableHeadings\[Rule]{Automatic,Automatic}]*)}], 
      Row[{"Curr node, newnode->", {ransamloop[[w]], newnode}}],
      Row[{"fuconodes->", fuconodes}],
      Row[{"activenodes->", activenodes}],
      Row[{"todel-> ", 
        "numtodel:" <> ToString[{todel1, todel2}] <> 
         " ", {If[loopnodecomplete, ransamloop[[w]], ""], 
         If[connodecomplete, newnode, ""]}}],
      Row[{"toadd->", toadd1, toadd2}],
      Row[{"samfwdcleaned->", samfwdcleaned}],
      Row[{"forward sample->", ransamloop[[w + 1 ;; iterations]]}]}]
    ]],
  (*loop print section in debug mode END  \
====================================== *)
  {w, iterations}](* end of main loop *); 
 Manipulate[
  Column[{Row[{"n=" <> ToString[n], "k=" <> ToString[k], 
      "\[Alpha]=" <> ToString[\[Alpha]value], Row[{"p=", pvalue}]}, 
     ", "], Row[{"MeanClusteringCoeff==" <> 
       ToString[
        Round[MeanClusteringCoefficient[AdjacencyGraph[matdyn[t]]], 
         0.01]], "L=MeanGraphDistance=" <> 
       ToString[
        Round[MeanGraphDistance[AdjacencyGraph[matdyn[t]]], 0.01]]}, 
     ", "], GraphPlot[matdyn[t], VertexLabeling -> True, 
     VertexCoordinateRules -> 
      Table[{s -> {Cos[\[Pi]/2 - (s - 1) (2 \[Pi])/Length[matdyn[t]]],
            Sin[\[Pi]/2 - (s - 1) (2 \[Pi])/Length[matdyn[t]]]}}, {s, 
         1, n + 1}][[All, 1]], ImageSize -> 600]}],
  {{t, iterations}, 0, iterations, 1, Appearance -> "Labeled"}, 
  TrackedSymbols :> {t}],
 Row[{Framed[
    Column[{Row[{"Graph construction computation"}, 
       Alignment -> Right, BaseStyle -> {Small, Italic}],
      Row[{Dynamic[
         ProgressIndicator[progr, {0, 100}, ImageSize -> {100, 7}]], 
        Pane[Dynamic[ToString[progr] <> "%"], ImageSize -> {30, 10}, 
         Alignment -> Center]}]}], FrameMargins -> {{5, 5}, {5, 5}}, 
    ContentPadding -> False, FrameStyle -> Red, Alignment -> Right],
   Control@{{\[Alpha]value, 0.75}, {0, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 
      4, 5, 8, 15}}, 
   Control@{{pvalue, 0.0001}, {10.^-9, 0.0001, 0.01, 0.1, 0.2, 0.3, 
      0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1}}, 
   Control@{{n, 20}, {10, 15, 20, 25, 30, 40, 50(* takes too long*)}},
    Control@{{k, Floor[n/4]}, {Max[2, Floor[n/10]], Floor[n/4], 
      Floor[n/2], Floor[(3 n)/4], n - 1}, ControlType -> Setter}}, 
  "   "], TrackedSymbols :> {\[Alpha]value, pvalue, n, k}]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.