Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to generate a version of Watts' model of cascading shocks to a network following the discussion on p. 20-25 of this article.

The basic premise is to use the RandomGraph function to generate a graph, then create a table of characteristics on from VertexList. The values in the table (x) are the vertex number, the 'persuadability' of the vertex, and the current state (0,1). Everyone is at 0 to start, and we put in an initial shock of 5 vertices going from 0 to 1. Then if the ratio of neighbours who are 1s is above the 'persuadability' number, that neighbour changes its state from 0 to 1.

UPDATE 2

The code below does what I want, in that it iterates the cascade and records each step, and also does it for different initial system shocks.

The code does most of what I originally asked about -

Ultimately I would like to be able to do the following:

  1. Use the code to store a table of x tables at each iteration. This could be used to animate the graph later (using HighlightGraph[SubGraph[...]]).

  2. Use the code to simulate various random shocks on the same graph by iterating until the table x does not change.

  3. Maybe a nice Manipulate interface at some point.

However, the evaluation takes about 20 minutes, and it also iterates a fixed number of times iter instead of sufficient for convergence. So I have some specific questions.

  1. Why is the assignment record[[l]]=x necessary in the Table loop?
  2. How can I get the Table to iterate till record[[l]]==record[[l-1]]? (Mike suggests FixedPointList, but I can't seem to get the syntax to work)
  3. How can I write the Do loop with the If statement so I can use ParallelDo or similar? (much like the answer here, but with conditions from two lists)
  4. What are some other ways to speed this up generally?

UPDATE 3

This part of the code takes about half a second (which seems slow, since it is called every iteration, or 750 times)

nextn = Table[
          Drop[VertexList[NeighborhoodGraph[y, neighbourVL[[i]]]], 1], {i,
            1, Length[neighbourVL]}];

This generates a table with each row a list of neighbours of the vertexes in the list neighbourVL. Is there a way to get this list with a faster method?

Thanks everyone.

  (* Note 20min evaluation time *)

y = RandomGraph[{100, 150}];
edges = EdgeList[y];
nodes = Partition[
   Flatten[Riffle[VertexList[y], Table[{RandomReal[0.7], 0}, {100}]]],
    3];
iter = 15;
record = Table[0, {iter}];

output = Table[

   shock = Transpose[{RandomSample[Range[100], 5], Table[3, {5}]}];
   x = ReplacePart[nodes, shock -> 1];

   Table[
    converted = Cases[x, {_, _, 1}];

    neighbourVL = 
     Flatten[Table[
       Drop[VertexList[NeighborhoodGraph[y, converted[[i, 1]], 1]], 
        1], {i, 1, Length[converted]}]];

    (* Vertex list of neighbours' first neighbours to be evaluated
with If function *)

    nextn = Table[
      Drop[VertexList[NeighborhoodGraph[y, neighbourVL[[i]]]], 1], {i,
        1, Length[neighbourVL]}];

    (* Number of neighbours of vertex 1's first neighbour who are 1s *)

        neighbourcount = 
     Table[Table[
       Count[x, {Take[nextn[[i, j]]], _, 1}], {j, 1, 
        Length[nextn[[i, All]]]}], {i, 1, Length[nextn]}];

    threshold = 
     N[Total[neighbourcount, {2}]/
       Table[Length[nextn[[i, All]]], {i, 1, Length[nextn]}]];

    Do[

     If[threshold[[i]] > x[[neighbourVL[[i]], 2]], 
       x = ReplacePart[x, {neighbourVL[[i]], 3} -> 1]];, {i, 
      Length[nextn]}];


    record[[l]] = x

    , {l, 1, iter}]

   , {50}];


Histogram[Total[output[[All, -1, All, -1]], {2}]]

Manipulate[
 Table[HighlightGraph[y, 
   Subgraph[y, Cases[output[[j, i]], {_, _, 1}][[All, 1]]]]
  , {j, 1, 6}], {i, 1, iter, 1}]
share|improve this question
    
kguler, for me the {Manipulate} seems to work with or without the {With}. Excuse the beginner question, but what is the purpose of defining l as itself? –  Cameron Murray Oct 23 '12 at 8:53
    
you are right, it does work with or without With. Sorry for the confusion. –  kguler Oct 23 '12 at 10:01
    
The link doesn't go to an article by Watt. I found this link (oz.stern.nyu.edu/phd04/watts.pdf) to a PNAS article, but I'd like to see yours. Important topic. –  George Wolfe Oct 23 '12 at 12:17
    
Thanks George. Your link is the original article. The link I have goes to a paper by Ormerod which simply had a very basic example of Watt's model. –  Cameron Murray Oct 24 '12 at 0:11
add comment

1 Answer 1

up vote 2 down vote accepted

I have code that works about 50x faster now. The trick was to use the GraphUtilities package function NeighborhoodVertices[] rather than VertexList[NeighborhoodGraph[]]. It sped up these two lines from about 0.5sec each to 0.007sec. Since these are called 1000 times or more, that's a big help.

I am still interested in a more functional way of doing this. I have two Tables nested with each other, with the inner Table running a Do loop each iteration. Any pointers to previous questions about converted code to functional and faster code is appreciated (I think I have read more of them). It is probably the inner Table function that could be a separate function to map over a table of shocked node lists (x).

It takes about 30sec now (40 shocks with 20 iterations each). Would be great to improve that.

Here is the code

shockNum = 40;
iter = 20;
y = RandomGraph[{100, 130}];
nodes = Partition[
   Flatten[Riffle[VertexList[y], Table[{RandomReal[0.6], 0}, {100}]]],
    3];

output = Table[
   shock = Transpose[{RandomSample[Range[100], 5], Table[3, {5}]}];
   x = ReplacePart[nodes, shock -> 1];

   Table[
    converted = Cases[x, {_, _, 1}];
    neighbourVL = 
     Flatten[Table[
       Drop[NeighborhoodVertices[y, converted[[i, 1]], 1], 1], {i, 1, 
        Length[converted]}]];
    nextn = 
     Table[Drop[NeighborhoodVertices[y, neighbourVL[[i]], 1]], {i, 1, 
       Length[neighbourVL]}];
    neighbourcount = 
     Table[Table[
       Count[x, {Take[nextn[[i, j]]], _, 1}], {j, 1, 
        Length[nextn[[i, All]]]}], {i, 1, Length[nextn]}];
    threshold = 
     N[Total[neighbourcount, {2}]/
       Table[Length[nextn[[i, All]]], {i, 1, Length[nextn]}]];
    Do[
     If[threshold[[i]] > x[[neighbourVL[[i]], 2]], 
       x = ReplacePart[x, {neighbourVL[[i]], 3} -> 1]];, {i, 
      Length[nextn]}];
    x
    , {iter}]
   , {shockNum}];

Manipulate[
 Column[{
   Histogram[Total[output[[All, i, All, -1]], {2}], 10, 
    PlotRange -> {{0, 100}, {0, shockNum/2}}, ImageSize -> 550, 
    AspectRatio -> 0.3, AxesOrigin -> {0, 0}],
   Partition[
     Table[HighlightGraph[y, 
       Subgraph[y, Cases[output[[j, i]], {_, _, 1}][[All, 1]]]]
      , {j, 1, 6}], 3] // Grid}]
 , {{i, 1, "Iteration"}, 1, iter, 1, Appearance -> "Labeled", 
  ContinuousAction -> False}, 
 FrameLabel -> {{None, None}, {None, 
    "Watt's model of cascading network failure"}}, 
 AppearanceElements -> All, LabelStyle -> Directive[Black, Medium]]

And here is the output.

enter image description here

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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