Table Of ContentRapid algorithm for identifying backbones in the two-dimensional percolation model
∗
Wei-Guo Yin
Department of Physics, Fudan University, Shanghai 200433, People’s Republic of China
and Physics Department, the Chinese University of Hong Kong, Hong Kong, People’s Republic of China
Ruibao Tao
1
CCAST (World Laboratory), P. O. Box 8730 , Beijing 100080, People’s Republic of China
0
and Department of Physics, Fudan University, Shanghai 200433, People’s Republic of China
0
2 (June 16, 1998; v2, September22, 1999; v3, February 1, 2008)
n Wepresentarapidalgorithmforidentifyingthecurrent-carryingbackboneinthepercolationmodel.
a It applies to general two-dimensional graphs with open boundary conditions. Complemented by
J the modified Hoshen-Kopelman cluster labeling algorithm, our algorithm identifies dangling parts
5 using their local properties. For planar graphs, it finds the backbone almost four times as fast as
1 Tarjan’sdepth-first-searchalgorithm,andusesthememoryofthesamesizeasthemodifiedHoshen-
Kopelman algorithm. Comparison with other algorithms for backboneidentification is addressed.
]
n
PACS numbers: 64.60.Ak, 02.70.-c, 05.10.-a
n
-
s The percolation model describes a system consisting algorithm [11] or the multigrid algorithm [14,15].
i of randomly distributed “conducting” cells and “isolat- The purpose of this paper is to present a fast algo-
d
. ing” cells [1]. When the density of the conducting cells rithm for d = 2 backbone identification. The backbone
t
a exceeds a threshold value pc, a cluster of connected con- is called the set of biconnected nodes in computer sci-
m ducting cells spans the system. This model has become ence. Forgeneralgraphs,heretoforethefastestalgorithm
- one of the most extensively studied statistical models to identify it is Tarjan’s depth-first-search algorithm [7]
d because of its simplicity and its often surprising appli- which runs in time O(N) for a graph of N nodes, pro-
n
cability. Its relevance to physics includes, for example, vided the number of edges meeting at each node is fi-
o
themetal-insulatortransitionobservedindisorderedsys- nite. For d = 2 planar graphs, our algorithm is almost
c
[ tems [2,3]. Numerical simulation is an important means four times as fast as Tarjan’s algorithm. We emphasize
in the solution of the model. For a complex system, here that our algorithm applies to general d = 2 graphs
2
which can be originated from the complexity of its ge- with open boundary conditions. There are other algo-
v
ometry or the mechanism underlying the random dis- rithms for backbone identification [5,6,8,10,12]. We will
9
6 tribution function, numerical simulation seems to be a first describe our algorithm and then compare it with
1 unique method available. To get rid off severe statisti- these known algorithms.
2 cal fluctuations brought by disorder in the percolation Below we give the definitions of all terminologies in-
1 model, numerical calculations, in general, need to be volved and illustrate them in a square lattice having L
0 done overa largenumber oflargescale ensembles. Find- siteshorizontallyandL+1sitesvertically(seeFigure1).
0
ing rapid algorithms is thus an interesting topic in this Betweenneighboringsites weinsert abond that is taken
/
t field [4,5,6,7,8,9,10,11,12,13,14,15,16]. to be open with an independent probability p, or closed
a
m Thecriticalbehaviorofpercolationiswellunderstood. with probability 1 − p. Obviously, here is defined the
Because of its relationship to the one-state Potts model percolationmodelforbonds[19]. Openlateralboundary
-
d [1], allcriticalexponents that havethermal analogueare conditionsareassumed,whereasallbondslocatingatthe
n exactly known for the dimension d = 2. Those having topandthebottomaresetopen. Apotentialdifferenceis
o no thermal analogue are usually estimated by numerical appliedbetweenthe topandthe bottom. Two neighbor-
c work. The most important of these are the backbone ing sites are connectedto eachotherif the bond between
:
v fractal dimension and the conductivity exponent. The them is open. The spanning cluster is the set of sites if
i backboneis defined asthe subsetofthe spanning cluster there is a path of open bonds from each of them to the
X
that carries current when a potential difference is ap- bottom and to the top. The spanning cluster consists of
r
plied between two sides of the system. Far away from the current-carrying backbone and non-current-carrying
a
p , the conductivity of the system can be accurately cal- danglingparts. Therearethreekindsofdanglingpartsof
c
culated within the effective-medium theory [2,17]; oth- the spanning cluster: dangling ends, dangling arcs, and
erwise, numerical studies are needed. The calculation danglingloops. Adangling endisconnectedtothe span-
of the conductivity on the current-carrying backbone is ning cluster by only one open bond (e.g., b in Figure 1);
considerably faster than on the spanning cluster [12,14] A dangling arc disconnects with the top (bottom) but is
becausethedensityofthebackboneinthespanningclus- connectedtothebottom(top)byatleasttwoopenbonds
ter substantially decreases as the system approaches to (e.g., c’s in Figure 1); A dangling loop is connectedby at
p [18]. Backbone identification is thus a necessary step least two open bonds to only one site of the spanning
c
beforeestimatingthe conductivity usingthe Lobb-Frank cluster (e.g., a in Figure 1). Here we name this site a
1
dangling-loop-connecting site,whichiscalledanarticula- the same label. This is the application of Corollary 1.
tion node in computer science. After all dangling parts Step 3foridentifyingdanglingarcs. Iftwoplaquettes
areremoved,theremainderofthespanningclusteristhe neighboringthetop(orthebottom)havethesamelabel,
backbone. then close all bonds between the two plaquettes. This is
OuralgorithmisbasedonJordan’s curve theorem: the application of Corollary 2.
Supposethereisaloopinatwo-dimensionalgraph,then Step 4 for identifying dangling loops. We look
theareainsidetheloopdisconnectswiththeareaoutside, for dangling-loop-connecting sites (a) using Corollary 3.
and vise versa. Thatis,aconnectsfouropenbonds,referredtoasbonds
To make use of the theorem,we introduce the concept 1,2,3,4, and a plaquette between bond 1 and 2 has the
of plaquette that is an as small as possible area enclosed same label as a plaquette between bond 3 and 4. Note
by a loop of bonds which can be either open or closed. that a’s themselves belong to the backbone. We use the
Forexample,anyoneofthesmallestsquareareasinFig- following trick to split dangling loops: We create a new
ure 1. We define that two nearest neighboring plaquettes site a′ and let bonds 2 and 3 be connected to a′ instead
are connected if the bond between them is closed. Note ofa. Asaresult,thesubgraphconnectedtobonds2and
that the definition of connectedness for plaquettes is op- 3 is disconnected to that connected to bonds 1 and 4.
posite to that for sites. Jordan’s curve theorem tells us Step 5. Computethe connectednessofsitesincluding
′
that any dangling part must be enclosed by a loop of those newly created sites (a). The resultant spanning
connectedplaquettesifitissplitfromthe spanningclus- cluster is exactly the backbone.
ter. Hence, by use of the connectedness of plaquettes, In this paper, the connectedness of plaquettes or sites
we are able to identify dangling parts according to the is calculated by using the modified Hoshen-Kopelman
followingthreecorollariesofthetheorem,whichdescribe cluster-labeling algorithm [4,14]. The data structure of
the local properties of the three kinds of dangling parts, this algorithm includes a bond status array that records
respectively. We will directly use them to design our al- thestatusofanybond[21],asitelabel arraythatrecords
gorithm. the label of any site (connected sites will have the same
Corollary 1. An open bond connects a dangling end label), and a bond-site array that records two end sites
to the spanning cluster if and only if two plaquettes be- of any bond. Our algorithm needs two extra arrays to
side the bond are connected. Forexample, b in Figure 1. record the information of plaquettes. One is a plaquette
Corollary 2. Iftwo non-neighboringplaquettes near- label array that records the label of any plaquette (con-
est neighboring the top (bottom) are connected, all sites nected plaquettes will have the same label), the other is
thatareconnectedto the top(bottom) bythe bonds be- a bond-plaquette array that records two side plaquettes
tween the two plaquettes belong to dangling arcs. For of any bond. Given this data structure, we can compute
example, c in Figure 1. boththe connectednessofsites andthe connectednessof
Corollary3. Letalldanglingendsandarcshavebeen plaquettes. Theprocessesofthealgorithmandtheusage
removed. Site a is a dangling-loop-connecting site if and ofthe arraysaresummarizedinFigure2. Onlythebond
only if (i) a connects at least four open bonds, referred status array is used in all steps. It is clear that if the
to as bonds 1,2,3,4, and (ii) among its neighboring pla- memory is the bottle-neck of calculations, the bond-site
quettes,somebetweentwoofthefourbonds(e.g.,bonds arrayandthebond-plaquettearraycanoccupythe same
1 and 2) is connected by a path of plaquettes to some memory,socanthesitelabelarrayandtheplaquettear-
others between the other two bonds (i.e., bonds 3 and ray. Thereasonsare(i)theyarenotusedsimultaneously
4). For example, a in Figure 1. a’s northern, eastern, and (ii) they are used sequentially. The bond-plaquette
southern, and western open bonds correspond to bonds arrayisusedinStep1-3whilethebond-sitearrayisused
1, 2, 3, 4, respectively. a’s northeast plaquette which is inStep4-5. TheplaquettelabelarrayisusedinStep1-4
betweenbonds 1and2 isconnectedby apathofplaque- while the site label array is used in Step 5.
ttes to a’s southwest plaquette which is between bonds Figure 3 shows one of the simplest d = 2 nonplanar
3 and 4. Note that here bonds 1 and 4 belong to the graphs: a square lattice with both nearest neighboring
backbone, bonds 2 and 3 belong to a dangling loop. (nn) bonds and next nearest neighboring (nnn) bonds.
Nowletusdescribetheprocedureofouralgorithm. In This type of nonplanar graphs is important because it is
accordance with open lateral boundary conditions, the often used to mimic the effects of higher dimensionality
left(right)plaquettesoftheleftmost(rightmost)column [15]. Hereaplaquetteisasmallesttriangulararea. While
of vertical bonds are set to be the same, respectively, as a nn bond still has two side plaquettes, a nnn bond has
showninFigure1. Thealgorithmiscarriedoutasfollows four which are divided into two pairs by another nnn
(see Figure 2). bond. If any pair of side plaquettes of a nnn bond have
Step 1. Compute the connectedness of plaquettes. the same label, the nnn bond can be closed in Step 2.
If the leftmost plaquette is connected to the rightmost The number of bonds (plaquettes) in this case is twice
plaquette by a path of plaquettes, the spanning cluster (four times) as large as that in the correspondingplanar
does not exist [20]. square lattice with nn bonds only.
Step 2 for identifying dangling ends. Sweep the bond TheefficiencyofouralgorithmisdemonstratedinFig-
graph to close any bond whose two-side plaquettes have ure 4 for square lattices with nn bonds only, where the
2
central-processing-unit(CPU) time neededfor the simu- most efficient algorithm computing the connectedness of
lation is plotted versus the number of sites. We consider sites, which is currently the modified Hoshen-Kopelman
bond percolation at p = p(b) = 0.5 [3] and site percola- algorithm. The modified Hoshen-Kopelman algorithm
c
tion [19] at p = p(s) = 0.592745(2) [22]. The results are was originally used to identify the spanning cluster [4].
c
Because of its simplicity, the geometrical properties of
obtained by averaging over 1000 configurations [21] for
the spanning cluster have been numerically studied on
linear size L = 100,200,300,500,1000,1500,2000,2500,
very large systems containing as much as 1011 sites [4].
respectively. Results for systems of larger size and more
Therefore,the inneroperationsofouralgorithmaresim-
configurations will be reported soon. The calculations
plerthanthoseofTarjan’sandhull-generating. Itshould
are carried out on a Pentium II/233 processor with Red
be made clear that the algorithms of Tarjan, matching,
HatLinuxrelease4.2operatingsystemandGNUproject
and burning can be used for arbitrary graphs; the hull-
F77Compiler(v0.5.18). AsshowninFigure4,thedotted
generating algorithm is valid for strictly planar graphs
lines are our fits to the data for three algorithms. Their
only; ours applies to general d=2 graphs with open lit-
slopes are 0.51, 1.21, and 4.62 for the modified Hoshen-
eral boundary conditions, and thus has a wider range of
Kopelman algorithm, our new algorithm, and Tarjan’s
application than the hull-generating algorithm.
depth-first-search algorithm, respectively. All the three
Theexactfindingofthespanningclusterandtheback-
algorithms run in time proportional to the system size.
bone enables us to calculate some scaling exponents for
Thetotaltimeneededtofindthebackboneisabouttwice
as many as that to compute the connected cluster, since percolation. The order parameter of percolation is P∞,
theratioofthesitesinthespanningclustertoallsitesin
weusethemodifiedHoshen-Kopelmanalgorithmtwicein
the formerandonceinthe latter. The presentalgorithm the system. Correspondingly, we define Q∞, the ratio of
thesitesinthebackbonetoallsitesinthesystem. These
finds the current-carryingbackbonealmostfour times as
quantities scale with linear size L as
fast as Tarjan’s depth-first-search algorithm. Note that
in our algorithm finding the backbone is not subject to P∞(L)∼L−β/ν, Q∞(L)∼L−q/ν, for L→∞ (1)
finding the spanning cluster first. An advantage of the
newalgorithmisthatitiseasytoprogramandmaintain where β, q, and ν are the critical exponents of P∞, Q∞,
because it involves only three local properties of dan- and the correlation length, respectively.
gling parts and uses the well-known modified Hoshen- We calculate these scaling exponents in our Monte
Kopelman algorithm. The codes of the algorithm are Carlo runs. Note that the following numerics are
available from us. listed for demonstration because we have not calcu-
There areother algorithmsforbackbone identification lated them on a world-record-breakingnumber of world-
[5,6,8,10,16]. A brief review of them was recently given record-breaking scale lattices. The largest size consid-
in Ref. [12]. Here we make a comparison between ours ered here is L = 2500 with 1000 configurations. Re-
andtheseknownalgorithm. The traditionallyusedalgo- sults for L = 4096 with 198470 configurations were
max
rithm by physicists is the burning algorithm [5,6], which reported in Ref. [12]. Calculations for systems of larger
is at least for large N much slower than Tarjan’s algo- size and more configurations are in progress. The data
rithm [12]. The matching algorithm with complexity are plotted in Figure 5 with very small standard errors
slight larger than N [O(N1.07) for d = 2] was used re- of the mean. The lines shown in the log-log plot are the
cently in literature [8]. For strictly planar graphs, the results of the weighted-least-squares fits, and yield the
hull-generating algorithm is even faster [9,10]. This al- value β/ν = 0.1068±0.0013 (close to the exact result
gorithmhas the same asymptotic complexityasTarjan’s β/ν =5/48[1,4]),andq/ν =0.3647±0.0039. Thecalcu-
algorithm, but it is roughly twice as fast and uses about latedexponentssatisfytheChayes’exponentinequalities
halfofthe memory,sinceitneedsonedatastructureless [18]:
and needs only one pass through all sites, instead of two
passesin Tarjan’s algorithm[12]. For planar graphs,our 2β ≤q ≤t, (2)
algorithm is almost four times as fast as Tarjan’s algo-
where t is the exponent of conductivity (the best esti-
rithm anduses the memoryofthe same size as the mod-
mate for two-dimensional systems known to us is t/ν =
ified Hoshen-Kopelman algorithm. For site percolation
0.9745±0.0015 [23]). The critical behavior of the back-
[19], a modified version of our algorithm uses the mem-
bone is consistently closer to that of the conductivity
oryofevensmallersizewherethebond-sitearrayandthe
bond-plaquette array are discarded [16]. Since the hull- thanthatofthespanningcluster. SinceQ∞/P∞vanishes
as p → p , the unknowns in the calculation of conduc-
generating algorithm, Tarjan’s, and ours have the same c
tivity onthe backbone is considerablyless thanthose on
asymptotic complexity,their difference in speed is deter-
the spanning cluster. As a result, the calculationof elec-
mined by the complexity of inner operations. Note that
tronicconductivityonthebackbonehasconsiderablyless
the identification of the backbone or the spanning clus-
statistical fluctuation and much reduced critical slowing
ter involves the global properties of sites. Our algorithm
down than on the spanning cluster [14].
hasanadvantagesuchthatitisabletoidentifydangling
Summarizing, we present a fast algorithm for identi-
parts using their local properties complemented by the
fying percolation backbones in general two-dimensional
3
graphs with open boundary conditions. Complemented [11] C.J.LobbandD.J.Frank,Phys.Rev.B30,4090(1984);
by the modified Hoshen-Kopelman algorithm, our algo- ibid. 37, 302 (1988).
rithm identifies dangling parts using their local proper- [12] P. Grassberger, Physica A 262, 251 (1999).
ties. Forplanargraphs,ouralgorithmisbyfarthefastest [13] M.E.J.NewmanandR.M.Ziff,Phys.Rev.Lett.85,4104
algorithm for backbone identification: it finds the back- (2000).
bone almost four times as fast as Tarjan’s depth-first- [14] R.G. Edwards, J.Goodman, andA.D.Sokal,Phys.Rev.
search algorithm, and uses the memory of the same size Lett.61,1333(1988);J.RugeandK.Stuben,inMultigrid
as the modified Hoshen-Kopelman algorithm. Methods,editedbyS.F.McCormick(SIAM,Philadelphia,
1987).
We thank Paul Bastiaansen for providing us with
[15] P.Bastiaansen andH.Knops,J.Phys.A30,1791(1997).
the code of Tarjan’s depth-first-search algorithm imple-
[16] Wei-GuoYin and RuibaoTao, Physica B 279, 84 (2000).
mentedbyJ.Goodman(Ref.[14])andM.E.J.Newman
[17] Wei-Guo Yin and Ruibao Tao, Phys. Rev. B 62, 550
for pointing out the original hull-generating algorithm
(2000).
(Ref. [9]) and his recent work (Ref. [13]). We are grate-
[18] J. T. Chayes and L. Chayes, Phys. Rev. Lett. 56, 1619
ful to our referees for their useful suggestions and point-
(1986); Comm. Math. Phys. 105, 133 (1986).
ing out recent important work done on the problem of
[19] Percolation,ofcourse,canbedefinedforsites.Inthiscase,
percolation backbones (Ref. [5,6,8,10,12]). W.G.Y. was
each site is taken to be occupied with probability p or va-
supportedbytheChinaPostdoctoralScienceFoundation cantwithprobability1−p,andabondisopenifandonly
and the Shanghai Postdoctoral Science Foundation, and
if both its end sites are occupied. But in the language of
is supported by CUHK 4288/00P 2160148. R.T. was graphtheorybondpercolationismoregeneralbecausethe
supported by the National Natural Science Foundation graphs created by site percolation are induced subgraphs
of China (No. 19774020),the foundation of Educational of the lattice.
Ministry of China (No. 99024610), and Shanghai Re- [20] It is not true for d = 2 percolation with periodic lateral
search Center of Applied Physics. boundary conditions.
[21] Configurations are produced by randomly producing the
bond status array: 1≡open and 0≡closed. Open bonds
appear with an independentprobability p.
[22] R.M. Ziff and B. Sapoval, J. Phys. A 19, L1169 (1986).
[23] J. M. Normand, H. J. Herrmann, and M. Hajjar, J. Stat.
Phys.52, 441 (1988).
∗
Permanent electronic address: [email protected].
[1] D. Stauffer and A. Aharony, Introduction to Percolation
Theory, 2nd edition (Taylor & Francis, Philadelphia, DA,
c
1994).
[2] S. Kirkpatrick,Rev.Mod. Phys. 45, 574 (1973).
[3] J. W. Essam, in: Phase Transitions and Critical Phenom-
a
ena,Vol.2,ed.C.DombandM.S.Green(Academic,New
York,1972) pp.218. 1
[4] J.HoshenandR.Kopelman,Phys.Rev.B14,3438(1976); b 4 2
H. Nakanishi and H. E. Stanley, ibid. 22, 2466 (1980); A. 3
Margolina,H.Nakanishi,D.Stauffer,andH.E.Stanley,J.
Phy.A17,1683(1984);D.C.Rapaport,J.Stat.Phys.66,
679(1992);H.GouldandJ.Tobochnik,AnIntroductionto
Computer Simulation Methods (Addison-Wesley,Reading,
MA, 1988), Part 2, pp.397-443; pp.637-645.
[5] M. D. Rintoul and H. Nakanishi, J. Phys. A 27, 5445
(1994); ibid. 25, L945 (1994).
[6] H. J. Herrmann, D. C. Hong, H. E. Stanley, J. Phys. A
17, L261 (1984); H. J. Herrmann and H. E. Stanley, ibid.
21, L829 (1988); S. Roux and A. Hansen, ibid. 20, L1281 FIG. 1. An illustration of percolation on a square lattice
(1987); M. Porto, A.Bunde, S.Havlin, and H.E. Roman, of linear size L = 10: site (filled circle), open bond (solid
Phys. Rev.E 56, 1667 (1997). line), closed bond (dotted line), plaquette (smallest square
[7] R.Tarjan, SIAMJ. Comput. 1, 146 (1972); A.V.Aho, J. area). aisadangling-loop-connecting site,bisanopenbond
E. Hopcroft, and J. D. Ullman, The Design and Analysis connecting a dangling end, and c are open bonds belonging
of Computer Algorithms (Addison-Wesley, Reading, MA, to a dangling arc. Dashed lines are some paths of connected
1974) pp.176-187. plaquettesenclosing some dangling parts.
[8] C. Moukarzel, Int.J. Mod. Phys. C 9, 887 (1998).
[9] R. M. Ziff, P. T. Cummings, and G. Stell, J. Phys. A 17,
3009 (1984).
[10] P. Grassberger, J. Phys. A 25, 5475 (1992).
4
30 Tarjan's depth-first-search
Processes Arrays DS) Hoshen-Kopelman
N new
O
STEP 1 C20
Plaquette Bond Status E
Connectivity S
OHHNSE| DaSngTlEinPg 2End Plaquette label U TIME (10
OKP Identification CP
PE STEP 3 0
MAL DIdaenngtilfiincga tAiornc Bond-Plaquette 0 2 4 6
N 10-6 LATTICE SITES
A
GL STEP 4
O Dangling Loop Band-Site
RI Identification
HT FIG. 4. CPU time of three algorithms for bond percola-
M STEP 5 tion at p = p(cb) = 0.5. The dotted lines are our fits to the
Site Site Label
Connectivity data. Their slopes are 0.51, 1.21, and 4.62 for the modified
Hoshen-Kopelman algorithm, the new algorithm, and Tar-
jan’s depth-first-searchalgorithm, respectively.
Backbone
FIG. 2. Flowchart of our algorithm. A line between an
array(parallelogram)andaprocess(rectangle)indicatesthat
thearrayisusedintheprocess. Ifthelineisdashed,thearray
changes in theprocess.
-0.5
ln(P)
-1.0
a b -1.5
-2.0
2
-2.5
1 3
ln(Q)
-3.0
4
5 6 7 8
ln(L)
d c
FIG. 5. Log-Log plot of Q∞ and P∞ for site percola-
FIG.3. A square lattice of linear size L=2 with crossing tion at p = p(s) = 0.5927. The dotted lines are the re-
c
bonds (ac and bd). Filled circles a,b,c,d denote sites. The sults of our weighted-least-squares fits against the functions
smallesttrianglularareas1,2,3,4denoteplaquettes. Bondac Q∞(L) = aL−q/ν and P∞ = bL−β/ν, giving values of
has two pairs of side plaquettes [(1,2) and (3,4)]. Bond bd q/ν =0.3647±0.0039 and β/ν =0.1068±0.0013.
hastwo pairsof sideplaquettes[(1,4) and(2,3)]. Ifanypair
ofsideplaquettesofacrossingbondhavethesamelabel, the
crossing bond can be closed.
5
