Finding a percolation path

Question

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]

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}]

Answer 1

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]

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]

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

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.

Answer 2

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}

Answer 3

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}]

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.

Answer 4

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 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]]}]