# Finding all minimal paths in site percolation?

I want to find all minimal site percolation paths from the left side to the right side, small demo grid.

A related question is here on finding bond percolation paths instead of minimal site percolation paths. A site is a neighbour to other site vertically or horizontally but not diagonally.

I want to use the FindShortestPath[], but I need to somehow get a GridGraph working: the basic idea is that unoccupied sites are removed but the vertex removal does not work as expected where I try to model the lattice as a grid graph, but I cannot understand why a vertex deletion results to a path graph, this puzzle moved here.

How can I find the minimal site percolation paths?

SeedRandom;
dimension = 100;
coDimension = 30;
percProbability = 0.7;
myData = Table[
Table[Boole[RandomReal[] < percProbability], {i, dimension}], {j,
coDimension}];
myData // MatrixPlot • Some people here would know how to convert this to a graph… I think using Sparsify?? You should make your question not an image though… Did you see mathematica.stackexchange.com/a/5156/1089 – chris Jul 9 '16 at 18:31
• @chris the image-based solution is computationally demanding on my comp, taking forever to compute :/ – hhh Jul 9 '16 at 18:36

Here's a solution using MorphologicalGraph[]:

SeedRandom;
dimension = 100;
coDimension = 30;
percProbability = 0.7;
myData = Table[Table[Boole[RandomReal[] < percProbability], {i, dimension}],
{j, coDimension}];
img = Binarize@Image@myData;


Now all you need to do is use FindShortestPath[]. For example, the shortest path from top-left to bottom-right corner:

g = MorphologicalGraph[img]
HighlightGraph[g, PathGraph[FindShortestPath[g, 1, Max@VertexList[g]]]] However, MorphologicalGraph[] includes corner neighbours, which we don't want in site percolation on a square lattice. A bit of digging turned up ImageMorphologicalOperationsDumpoMorphologicalGraph[] as the function behind MorphologicalGraph. By adapting this function to ignore the thinning operation and only use corner neighbours, you can get the appropriate graph. The adapted function is at the bottom of the post.

First, let's get the left and right vertices out from your data:

getIndex[sites_, site_] := Position[sites, _?(# == site &)]
getLeftAndRightVertices[data_] :=
Module[{sites, leftsites, rightsites},
sites = Position[data, _?(# == 1 &)];
leftsites = Select[sites, #[] == 1 &];
rightsites = Select[sites, #[] == Last@Dimensions@data &];
{Flatten[getIndex[sites, #] & /@ leftsites],
Flatten[getIndex[sites, #] & /@ rightsites]}]

{leftvertices, rightvertices} = getLeftAndRightVertices[myData];


Now you can find the shortest path between any left vertex and any right vertex. To get the overall shortest path from left to right, you can do the following, with a warning when no path can be found.

g2 = myMorphologicalGraph[img, VertexCoordinates -> Automatic];

allpairs = Tuples[{leftvertices, rightvertices}];
allpaths = Quiet@MapThread[FindShortestPath[g2, #1, #2] &, Transpose@allpairs];
pathlengths = Length@# & /@ allpaths;

(* Workaround to deal with unconnected components which *)
(* give lengths of 0. The minimum possible path length *)
(* is of course == dimension *)
If[Max@pathlengths >= dimension,
sortlengths = Ordering[pathlengths];
pos = FirstPosition[pathlengths[[sortlengths]], _?(# >= dimension &)];
shortestpair = Flatten@allpairs[[sortlengths]][[pos]];
shortestpath = Flatten@allpaths[[sortlengths]][[pos]];
Show[img, HighlightGraph[g2, PathGraph[shortestpath]]],
(* Print warning *)
Print["No connected path found"]]

(* Left = 770, Right = 1242 *)
(* Length = 128 *) Code for myMorphologicalGraph[]:

myMorphologicalGraph[skeleton_,
opts : OptionsPattern[MorphologicalGraph]] :=
Module[{vertices, vertexComponents, vertexCoordinates,
loopEdges, cleanEdges, extraEdges, allEdges,
(* Use the direct binarized image *)
vertices = skeleton;
vertexComponents =
Replace[ImageData[vertices], 1 :> ++vertexCount, {2}];
vertexCoordinates = OptionValue[VertexCoordinates];
If[vertexCoordinates === Automatic,
vertexCoordinates =
ComponentMeasurements[vertexComponents, "Centroid"]];
ImageMorphologicalOperationsDumpConstrainedMComponents[
ImageSubtract[skeleton, vertices],
Dilation[vertices, CrossMatrix]];
Replace[linkComponents, Except[0, n_] :> n + vertexCount, {2}];
all4all =
"Neighbors", CornerNeighbors -> False];
vertex4all = Select[all4all, First[#1] <= vertexCount &];
vertex4vertex =
DeleteCases[vertex4all, _?(#1 > vertexCount &), {3}];
vertex4link = DeleteCases[vertex4all, _?(#1 <= vertexCount &), {3}];
(* Set corner neighbours to False here *)
all8all =
"Neighbors", CornerNeighbors -> False];
vertex8all = Select[all8all, First[#1] <= vertexCount &];
vertex8vertex =
DeleteCases[vertex8all, _?(#1 > vertexCount &), {3}];
vertex8link = DeleteCases[vertex8all, _?(#1 <= vertexCount &), {3}];
link8all = Select[all8all, First[#1] > vertexCount &];
redundantEdges =
ImageMorphologicalOperationsDumpsortEdges[
DeleteCases[
ImageMorphologicalOperationsDumpgrowEdges[
ImageMorphologicalOperationsDumpgrowEdges[
ImageMorphologicalOperationsDumptoEdges[vertex4vertex],
edge[_, v_, _, v_] | edge[v_, _, _, v_]][[All, {1, 3, 4}]]];
directEdges =
Select[ImageMorphologicalOperationsDumptoEdges[vertex8vertex],
OrderedQ];
ImageMorphologicalOperationsDumpgrowEdges[
ImageMorphologicalOperationsDumptoEdges[vertex8link],
loopEdges =
Cases[Tally[linkedEdges, #1[] === #2[] &], {e_, 1} -> e];
loopEdges =
Pick[loopEdges,
cleanEdges =
extraEdges =
ImageMorphologicalOperationsDumpgrowEdges[
ImageMorphologicalOperationsDumpgrowEdges[link8vertex,
Select[Cases[
ImageMorphologicalOperationsDumpgrowEdges[
extraEdges =
ImageMorphologicalOperationsDumpsortEdges[
DeleteCases[extraEdges,
Alternatives @@
Cases[extraEdges,
edge[v_, l1_, l2_, v_] -> edge[_, l1, l2, _]]]];
allEdges =
Apply[UndirectedEdge,
directEdges \[Union] cleanEdges[[All, {1, -1}]] \[Union]
extraEdges[[All, {1, -1}]], {1}];
edgeWeights = OptionValue[EdgeWeight];
If[edgeWeights === Automatic,
edgeWeights =
Apply[UndirectedEdge, cleanEdges[[All, {1, -1}]], {1}] ->
1 + (cleanEdges[[All, 2]] /. linkWeights)]],
Thread[Apply[UndirectedEdge, extraEdges[[All, {1, -1}]], {1}] ->
3]]; edgeWeights =
Replace[allEdges, Dispatch[edgeWeights], {1}]];
Graph[allEdges, VertexCoordinates -> vertexCoordinates,
EdgeWeight -> edgeWeights,
Sequence @@
FilterRules[opts,
DeleteCases[
Options[MorphologicalGraph], (VertexCoordinates -> _) | \
(EdgeWeight -> _)]]]];

• Did you try to solve the shortest path between all left vertices and all right vertices? For me, the computer freezes in trying to do it. – hhh Jul 9 '16 at 21:20
• @hhh did you try the update to see if your computer can solve it in a suitable amount of time? It's pretty quick on my PC. – dr.blochwave Jul 10 '16 at 11:07
• Yes it works, nice one +1! – hhh Jul 10 '16 at 11:12
• Hmmm, I tried this with different lattice sizes and for some reason, it is unable to find the minimal path. Have you tried this with different lattice sizes? – hhh Jul 10 '16 at 12:04
• @hhh issue fixed I think, now issues a warning if no path found. – dr.blochwave Jul 10 '16 at 12:19

Solution based on the GridGraph

SeedRandom;
dimension = 20;
coDimension = 30;
percProbability = 0.7;
deleteMe =
Pick[Table[i, {i, dimension*coDimension}],
Table[RandomReal[] > percProbability, {i,
dimension*coDimension}]];
G = GridGraph[{dimension, coDimension}, VertexLabels -> "Name",
G = SetProperty[G, VertexCoordinates -> GraphEmbedding[G]];
H = VertexDelete[G, deleteMe]
FindShortestPath[H, 1, 600]
HighlightGraph[H, PathGraph[%]] that finds the shortest path from the site 1 to the site 600. And next I want to find the shortest path from the left side to the right side

rightSide = Complement[Table[i, {i, 581, 600}], deleteMe]
Table[FindShortestPath[H, 1, i], {i, rightSide}]
Table[Length[FindShortestPath[H, 1, i]], {i, rightSide}]
shortest = Table[FindShortestPath[H, 1, i], {i, rightSide}][];
HighlightGraph[H, PathGraph[shortest]] that is the shortest path of length 47 between the vertex 1 and the left side. Next I need to do this over each site on the right side to find the shortest path between the right side and the left side

paths = Table[
Table[FindShortestPath[H, j, i], {i, rightSide}], {j, leftSide}];
pathLengths = Table[
Table[Length[paths[[k]][[h]]], {k, Length[paths]}],
{h, Length[paths[[]][]]}]
Histogram[pathLengths, 50]
pathLengths // Max
FindShortestPath[H, rightSide[], leftSide[]]
HighlightGraph[H, PathGraph[%]]
Pick[pathLengths, pathLengths // Positive] where the zero entries tell me that there are sites from which no path to the other side. So we need to look for positive entries only to find the minimal paths and we found it to be of 33 size. The largest pairwise minimal path length is 53 and it turns out to be between 1-596.

FindShortestPath[H, rightSide[], leftSide[]]
HighlightGraph[H, PathGraph[%]] 