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 would like to examine percolation on a random lattice. To be exact, I wish to find the minimum length of a 'bond' needed such that the leftmost site can be connected to the rightmost site.

Here is an example of the lattice:

randPts = Table[RandomReal[{-10, 10}, 2], {200}]; 
randPlot = ListPlot[randPts, 
                PlotStyle -> {PointSize[0.0125]}, 
                PlotRange -> {{-10, 10}, {-10, 10}}, 
                AspectRatio -> 1, 
                Frame -> True]

random lattice

I have tried for a while to get this but have not had success. The basic plan was:

  1. Define a bond length $R$

  2. Look at each site one at a time. If another site(s) is within $R$ of a site, they will be in the same cluster. Each site will be in a cluster of 1 or more (obviously the larger $R$ chosen, the larger each cluster size)

  3. Take a site. Does it bond with other sites? If so then combine the two clusters together.

  4. Repeat step 3 for all sites.

  5. At the end ask if the leftmost cite and the rightmost sites are included in the conglomerate cluster. If so, percolation has occurred.

  6. Decrease $R$ and start over again until a threshold is found.

I think I am stuck somewhere in the step 3,4 area. Here is some of what I've tried: I have defined a module to find the distance between a site, j, and its nearest neighbor. The table, t, gives distance between j and all other sites:

minD[j_] := 
  Module[{},
    t = Table[{randPts[[i]], 
              Sqrt[(randPts[[j, 1]] - randPts[[i, 1]])^2 + (randPts[[j, 2]] - 
                 randPts[[i, 2]])^2]}, 
             {i, 1, Length[randPts]}];

    For[i = 1, i < Length[t] + 1, i++, 
      If[t[[i, 2]] == RankedMin[t[[All, 2]], 2], 
        coord[j] = t[[i, 1]] ]];
    Return[{coord[j]}];
  ];

This module takes the table of distances and picks out ones that are within the chosen bonding radius (1.5 here. the y>0 condition to so to not count the same site):

  cluster[k_] := 
    Module[{},
      minD[k];
      Return[
        Table[Cases[t, {x_, y_} /; y < 1.5 && y > 0][[i]][[1]], 
              {i, 1, Length[Cases[t, {x_, y_} /; y < 1.5 && y > 0]]}]];
    ]

So cluster[k] gives the sites within the cluster that is centered at site k. Now combining these clusters is what I am having a problem with. My idea was to start with a site and its cluster; find out what clusters that cluster intersects with and continue. I was not able to implement this correctly.

Another way to visualize or maybe solve the problem is in terms of increasing the site radius at each site until a percolation network is achieved:

 randMovie = 
   Manipulate[
    ListPlot[randPts, 
        PlotStyle -> {PointSize[x]}, 
        PlotRange -> {{-10, 10}, {-10, 10}}, AspectRatio -> 1, 
        Frame -> True], 
    {x, 0.00, 0.12, 0.002}]

enter image description here

share|improve this question
    
Incidentally, randPts can be written RandomReal[{-10, 10}, {200, 2}] –  Mr.Wizard May 4 '12 at 22:53
    
That part is now copyable –  BeauGeste May 5 '12 at 3:51
2  
Nice question. Please note that the purpose of Return is not really returning values from a Module. Instead of using Module[{}, ...; Return[result]; ] simply use Module[{}, ...; result] (not that lack of the final semicolon). In general, you almost never need to use Return in Mathematica, and when you do, it's good to be aware of some unusual behaviour ... –  Szabolcs May 5 '12 at 9:49

4 Answers 4

up vote 33 down vote accepted

A percolation network is just a kind of network, so I went in the direction of proposing a graph-theoretic approach. You seem to be measuring distances between nodes multiple times, but given the points don't move, you need only do it once:

ed = Outer[EuclideanDistance, randPts, randPts, 1];

You can get the positions of the nodes you are trying to connect like so:

leftmost = Position[randPts, {Min[randPts[[All, 1]] ], _}][[1, 1]]

rightmost = Position[randPts, {Max[randPts[[All, 1]] ], _}][[1, 1]]

Here is an auxiliary function that determines which nodes are no more than r distance from each other. I exclude zero distances to avoid the complication of self-loops.

linked[mat_?MatrixQ, r_?Positive] := Map[Boole[0 < # < r] &, mat, {2}]

It is easy to use this auxiliary function to create an adjacency matrix which can be visualised with the correct coordinates using the VertexCoordinates option.

gg = AdjacencyGraph[linked[ed, 2.], VertexCoordinates -> randPts]

Graph visualisation

Finding out whether the left-most and right-most points are connected is a matter of determining if FindShortestPath yields a non-empty result.

FindShortestPath[gg, leftmost, rightmost]
(* ==> {56, 16, 126, 156, 142, 174, 65, 49, 23, 88, 6, 45, 122, 68, 131, 139, 80} *)

Let's put all this together. I am going to build the option to test if the network is a percolation network in the same function that visualises the network.

Options[isPercolationNetwork] = {ShowGraph -> False}

isPercolationNetwork[points : {{_?NumericQ, _?NumericQ} ..}, 
  r_?Positive, opts : OptionsPattern[]] :=
  Module[{ed = Outer[EuclideanDistance, points, points, 1], 
   leftmost =  Position[points, {Min[points[[All, 1]] ], _}][[1, 1]], 
   rightmost = Position[points, {Max[points[[All, 1]] ], _}][[1, 1]]},
  With[{gg =  AdjacencyGraph[linked[ed, r], VertexCoordinates -> points]},
   If[OptionValue[ShowGraph],
    HighlightGraph[gg, PathGraph[FindShortestPath[gg, leftmost, rightmost]]], 
    Length[FindShortestPath[gg, leftmost, rightmost] ] > 1]]
  ]

If the option ShowGraph is True, it shows the graph and the connecting path; if it is False, it just returns True or False.

isPercolationNetwork[randPts, 2., ShowGraph -> True]

enter image description here

It is pretty straightforward to put all this together to find the minimum distance to create a percolation network.

minimumPercolationNetwork[points:{{_?NumericQ, _?NumericQ}..}, r0_?Positive] :=
 Module[{r = r0},
  While[isPercolationNetwork[randPts, r], r = r - 0.01]; 
  Print[r + 0.01]; 
  isPercolationNetwork[points, r + 0.01, ShowGraph -> True] ]

And the result:

minimumPercolationNetwork[randPts, 3.]
1.97

enter image description here

Execution is reasonably fast: Timing of the above example was a bit above 6s on my machine, but it depends on the initial value you pick for r.

share|improve this answer
2  
@BeauGeste, by the way, thank you for asking this question, I quite enjoyed working out an answer. Even after so many years using Mathematica, I'm amazed that so much is possible using so little code. I am not even a graph theorist and have never studied the topic. –  Verbeia May 5 '12 at 5:36
    
excellent answer. This was along the lines of what I was attempting though I was not aware of all the graph theory functions until now (I was using Ifs and For statements, ugly I know). It really is amazing how simple it can be done (for you experts :)) in Mathematica compared to coding it –  BeauGeste May 5 '12 at 18:14
2  
You can also do binary search on the r, reducing the complexity a bit more. –  Paxinum May 9 '12 at 20:33
1  
There is a really great example of Bond Percolation in Ch. 1 of Paul Wellin's book Programing with Mathematica. The example in the book uses GridGraph and RandomVariate[] with a BernoulliDistribution[] to assign bonds. It's quite clever. –  brown.2179 Jan 6 at 16:39

As it might be interest to others than me, it seems a generalization to 3D of @Verbeia's post would be

linked[mat_?MatrixQ, r_?Positive] := Map[Boole[0 < # < r] &, mat, {2}]
Options[isPercolationNetwork] = Flatten[{ShowGraph -> False, Options[HighlightGraph]}];

isPercolationNetwork[points : {{_?NumericQ, _?NumericQ, _?NumericQ} ..}, r_?Positive, 
opts : OptionsPattern[]]:= 
Module[{ed = Outer[EuclideanDistance, points, points, 1],
leftmost = Position[points, {Min[points[[All, 1]]], _, _}][[1, 1]],
rightmost = Position[points, {Max[points[[All, 1]]], _, _}][[1, 1]]}, 
With[{gg = AdjacencyGraph[linked[ed, r], 
   VertexCoordinates -> points /. {_, y_, z_} -> {y, z}]}, 
If[OptionValue[ShowGraph], 
 HighlightGraph[gg,PathGraph[FindShortestPath[gg, leftmost, rightmost]], 
  Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]], 
Length[FindShortestPath[gg, leftmost, rightmost]] > 1]]];

minimumPercolationNetwork[
points : {{_?NumericQ,_?NumericQ,_?NumericQ}..},r0_?Positive,opts: OptionsPattern[]]:= 
Module[{r = r0}, 
While[isPercolationNetwork[points, r],r =r-0.01];Print[r + 0.01];
isPercolationNetwork[points, r + 0.01, ShowGraph -> True, 
Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]]]

so that

 randPts = RandomReal[{0, 1}, {250, 3}];

and

 minimumPercolationNetwork[randPts, 0.2]

produces

 0.16

a 3D percolating graph

A version of the code which deals with different percolation directions, takes graph options and works in 2 and 3D is given below

linked[mat_?MatrixQ, r_?Positive] :=  Map[Boole[0 < # < r] &, mat, {2}]
Options[isPercolationNetwork] = 
Flatten[{ShowGraph -> False, PercolationDirection -> 1,  
Options[HighlightGraph]}];
isPercolationNetwork[points : {{_?NumericQ, _?NumericQ} ..}, 
r_?Positive, opts : OptionsPattern[]] := 
Module[{ed = Outer[EuclideanDistance, points, points, 1], leftmost, 
rightmost},
If[OptionValue[PercolationDirection] == 1,
leftmost = Position[points, {Min[points[[All, 1]]], _}][[1, 1]];
rightmost = 
 Position[points, {Max[points[[All, 1]]], _}][[1, 1]];,
leftmost = Position[points, {_, Min[points[[All, 2]]]}][[1, 1]];
rightmost = Position[points, {_, Max[points[[All, 2]]]}][[1, 1]];
];
With[{gg = 
  AdjacencyGraph[linked[ed, r], VertexCoordinates -> points]}, 
If[OptionValue[ShowGraph], 
 HighlightGraph[gg, 
  PathGraph[FindShortestPath[gg, leftmost, rightmost]], 
  Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]], 
 Length[FindShortestPath[gg, leftmost, rightmost]] > 1]]];

isPercolationNetwork[
points : {{_?NumericQ, _?NumericQ, _?NumericQ} ..}, r_?Positive, 
opts : OptionsPattern[]] := 
Module[{ed = Outer[EuclideanDistance, points, points, 1], leftmost, 
rightmost},
Which[OptionValue[PercolationDirection] == 1,
leftmost = Position[points, {Min[points[[All, 1]]], _, _}][[1, 1]];
rightmost = 
 Position[points, {Max[points[[All, 1]]], _, _}][[1, 1]];,
OptionValue[PercolationDirection] == 2,
leftmost = Position[points, {_, Min[points[[All, 2]]], _}][[1, 1]];
rightmost = 
 Position[points, {_, Max[points[[All, 2]]], _}][[1, 1]];,
OptionValue[PercolationDirection] == 3,
leftmost = Position[points, {_, _, Min[points[[All, 3]]]}][[1, 1]];
rightmost = 
 Position[points, {_, _, Max[points[[All, 3]]]}][[1, 1]];
];
With[{gg = 
  AdjacencyGraph[linked[ed, r], 
   VertexCoordinates -> points /. {x_, y_, z_Real} -> {x, y}]},
If[OptionValue[ShowGraph],
 HighlightGraph[gg, 
  PathGraph[FindShortestPath[gg, leftmost, rightmost]]
  (*GraphPlot3D[ggh,VertexCoordinateRules-> 
  Thread[Range[Length[points]]->points],Axes->True,AxesLabel->{x,
  y,z},ViewPoint->{0,0,500}]*)
  , 
  Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]], 
 Length[FindShortestPath[gg, leftmost, rightmost]] > 1]]];
Clear[minimumPercolationNetwork];
Options[minimumPercolationNetwork] = 
Flatten[{ShowGraph -> True, PercolationDirection -> 1,  
Options[HighlightGraph]}];
minimumPercolationNetwork[points : {{_?NumericQ, _?NumericQ} ..}, 
r0_?Positive, opts : OptionsPattern[]] :=
Module[{r = r0},
While[isPercolationNetwork[points, r, 
PercolationDirection -> OptionValue[PercolationDirection]], 
r = r - 0.01];
{r + 0.01,
isPercolationNetwork[points, r + 0.01, ShowGraph -> True, 
 PercolationDirection -> OptionValue[PercolationDirection], 
 Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]] // 
Rasterize[#, ImageResolution -> 150] &}]

 minimumPercolationNetwork[
 points : {{_?NumericQ, _?NumericQ, _?NumericQ} ..}, r0_?Positive, 
 opts : OptionsPattern[]] :=
 Module[{r = r0},
 While[isPercolationNetwork[points, r, 
 PercolationDirection -> OptionValue[PercolationDirection]], 
 r = r - 0.01];
 {r + 0.01,
 isPercolationNetwork[points, r + 0.01, ShowGraph -> True, 
 PercolationDirection -> OptionValue[PercolationDirection],
 Sequence @@ FilterRules[{opts}, Options[HighlightGraph]]] // 
 Rasterize[#, ImageResolution -> 150] &}] 

 randPts = RandomReal[{0, 1}, {150, 3}];

 Column[{minimumPercolationNetwork[randPts, 
 1.5/(Length[randPts])^(1/3), Frame -> True, 
 PercolationDirection -> 1][[2]],
 minimumPercolationNetwork[randPts, 1.5/(Length[randPts])^(1/3), 
 Frame -> True, PercolationDirection -> 2][[2]],
 minimumPercolationNetwork[randPts, 1.5/(Length[randPts])^(1/3), 
 Frame -> True, PercolationDirection -> 3][[2]]}]
share|improve this answer

An image-based method ... just a curiosity:

r = 10; (*half range*)
i = step = 1/100;
rndpts = RandomReal[{-r, r}, {200, 2}];

l = Graphics[{Thickness[.001 r], Line@{{{-r, -r}, {r, -r}}, {{r, r}, {-r, r}}}}];
lPlot[i_] := ListPlot[rndpts, PlotStyle -> {Black, PointSize[i/(2 r)]}, 
                              PlotRange -> {{-r, r}, {-r, r}}, 
                              AspectRatio -> 1, Axes -> False];
t[i_] := MorphologicalComponents[ColorNegate@Binarize@Rasterize@Show[lPlot[i], l]];

(* Now loop until the image top and bottom rows are connected *)
While[(mem = t[i])[[1, 1]] != mem[[-1, 1]], i += (r step)];
{i, t[i] // Colorize}

enter image description here

share|improve this answer
    
That's pretty neat. Can you elaborate on your While statement? What is that doing to specify the percolation criterion? –  BeauGeste May 5 '12 at 4:45
    
@BeauGeste It is simply testing that the first pixel of the first row is in the same morphological component than the first pixel of the last row, so there is a connection between the top and bottom. Please note that I draw a line at the top and bottom of the graph to be able to test connectivity. –  belisarius May 5 '12 at 4:59
    
+1, I was thinking of doing the same when I read the question –  Szabolcs May 5 '12 at 9:40

I learned about this technique from Fred Simons on MathGroup, in a thread about computing connected components in graphs. You'll find the full discussion thread here.

Let's first create the sample dataset:

pts = RandomReal[10 {-1, 1}, {200, 2}];

ListPlot[pts, AspectRatio -> Automatic, 
 Epilog -> {Red, Point[pts[[63]]], Point[pts[[90]]]}]

Then let's compute a distance matrix between points:

dst = Outer[EuclideanDistance, pts, pts, 1]; // Timing

(If you wish, you can speed this up by not computing every distance twice. I chose to keep the code simple.)

Like @Verbeia, I chose to use a graph-apporach. Let's create the set of possible edges in the graph and sort them by length.

edges = Subsets[Range@Length[pts], {2}];
edges = SortBy[edges, Extract[dst, #] &];

Let' choose the leftmost and rightmost points and name their indices start and end:

start = First@Ordering[pts[[All, 1]], 1];
end = First@Ordering[pts[[All, 1]], -1];

And now use Fred's solution with a little modification:

idx = Module[{f}, 
       Do[
        Set @@ f /@ (edges[[i]]); 
        If[f[start] === f[end], Return[i]], 
        {i, Length[edges]}]]

idx will give the edge of length $R$ (i.e. the minimal length edge that needs to be included). In my case this length was 2.27:

Extract[dst, edges[[idx]]]

(* ==> 2.27273 *)

Here's a Manipulate that'll keep adding edges one by one, in order or increasing length, until we reach percolation. The leftmost and rightmost vertices are highlighted in red.

Manipulate[
 HighlightGraph[
  Graph[Range@Length@pts, UndirectedEdge @@@ Take[edges, i], 
   VertexCoordinates -> pts], {start, end}], {i, 1, idx, 1}]

Mathematica graphics

If the performance of this solution is not good enough, you can speed it up a little bit using the method I described in this MathGroup post. The total running time for 200 points is ~0.2 seconds on my (slow) computer.

share|improve this answer
    
+1 I was writing a similar solution to this but you beat me. Serves me right for having lunch. –  Heike May 5 '12 at 11:12
    
@Heike I almost went away to do something else for a while but I had a feeling that I needed to finish fast :-) –  Szabolcs May 5 '12 at 11:18

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.