(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "NeXT Mathematica Notebook Front End Version 2.2"; NeXTStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L1, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L1, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L1, a15, 14, "Courier"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L1, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L1, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L1, 12, "Courier"; ; fontset = name, inactive, noPageBreakInGroup, nohscroll, preserveAspect, M7, italic, B65535, L1, 10, "Times"; ; fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 12, "Times"; ; fontset = leftheader, 12; fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic, L1, 12, "Times"; ; fontset = leftfooter, 12; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Courier"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; paletteColors = 128; showRuler; automaticGrouping; currentKernel; ] :[font = title; inactive; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] Simulating Experiences: Excursions in Programming ;[s] 2:0,0;1,1;50,-1; 2:1,12,9,Times,1,14,0,0,0;1,21,16,Times,1,24,0,0,0; :[font = subtitle; inactive; preserveAspect; fontSize = 14; fontName = "Times"] ``Forest Fires and Reforestation'' :[font = subtitle; inactive; preserveAspect] Mathematica in Education Vol.3 No.2 Spring 1994 (c) TELOS/Springer-Verlag ;[s] 2:0,0;24,1;76,-1; 2:1,17,13,Times,3,18,0,0,0;1,16,12,Times,1,18,0,0,0; :[font = subsubtitle; inactive; preserveAspect] by Richard J. Gaylord Department of Materials Science University of Illinois at Urbana-Champaign gaylord@ux1.cso.uiuc.edu and Kazume Nishidate Department of Electrical and Electronic Engineering Iwate University, Faculty of Engineering Morioka 020 JAPAN nisidate@wriron1.s.kanazawa-u.ac.jp :[font = subsubtitle; inactive; preserveAspect] edited by Richard J. Gaylord :[font = section; inactive; Cclosed; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] Introduction :[font = text; inactive; preserveAspect; endGroup] Turbulent cascading has been suggested as the underlying cause of a wide variety of phenomena, including geological upheavals (eg., volcanic eruptions and earthquakes), species extinction during evolution, and fluid turbulence. A simple one-dimensional probabilistic cellular automaton, known as the forest fire model, displays this behavior. :[font = section; inactive; Cclosed; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] The Forest Fire CA :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] System :[font = text; inactive; preserveAspect] The forest fire CA employs a one-dimensional lattice of length n, with periodic boundary conditions. Lattice sites may have a value of 0, 1 or 2, where 0 represents an empty site (or hole), 1 represents a healthy tree (or tree), and 2 represents a burning tree. A forest is a set of contiguous sites (ie., a connected segment) with value 1 or 2. A forest preserve consists of forests separated by gaps (or deserts) consisting of clusters of connected holes. The system evolves in a specified number of time steps, in each of which entire forests of trees can catch fire and burn down and trees can sprout on individual empty sites. :[font = text; inactive; preserveAspect; endGroup] (Note: having all of the trees in a forest burn down in a single time step while trees grows independently of one another provides a separation between the time scales for the processes of deforestation and reforestation). :[font = subsection; inactive; preserveAspect; startGroup] Algorithm :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Statement and Implementation :[font = text; inactive; preserveAspect] (1) A forest preserve of length n, consisting of emtpy sites and tree sites is created using :[font = input; preserveAspect] forestPreserve = Table[Random[Integer], {n}] :[font = text; inactive; preserveAspect] All of the sites in the forest preserve are updated in each time step, according to the following sequence of steps: :[font = text; inactive; preserveAspect] (2a) Trees catch fire with probability f and empty sites sprout trees with probability, p, using a set of transformation rules to turn 0's into 1's with probability p, and 1's into 2's with probability, f :[font = input; preserveAspect] treeGrowIgnite = forestPreserve /. {0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]} :[font = text; inactive; preserveAspect] (2b) All tree sites adjacent to an ignited tree site (ie., trees in the same forest ) ignite. This accomplished using the repeated application of a set of transformation rules :[font = input; preserveAspect] forestIgnite = treeGrowIgnite //. {{a___, 2, 1, b___} -> {a, 2, 2, b}, {a___, 1, 2, b___} -> {a, 2, 2, b}, {2, c___, 1} -> {2, c, 2}, {1, c___, 2} -> {2, c, 2}} :[font = text; inactive; preserveAspect] (note: The first two transformation rules change sequences in the list having the form .., 1, 2,.. and .., 2, 1, .. to .., 2, 2, .. . and the last two transformation rules implement periodic boundary conditions by treating trees at each end of the system as belonging to the same forest. :[font = text; inactive; preserveAspect] (2c) All ignited trees burn down. :[font = text; inactive; preserveAspect] forestNew = forestIgnite /. 2 -> 0 :[font = text; inactive; preserveAspect] The sequence of steps 2a-c can be combined in an anonymous function which can be applied to any forest preserve configuration :[font = input; preserveAspect] pyro = (treeGrowIgnite = # /. {0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]}; forestIgnite = treeGrowIgnite //. {{2, c___, 1} -> {2, c, 2}, {1, c___, 2} -> {2, c, 2}, {a___, 2, 1, b___} -> {a, 2, 2, b}, {a___, 1, 2, b___} -> {a, 2, 2, b}}; forestNew = forestIgnite /. 2 -> 0)& :[font = text; inactive; preserveAspect] where # represents the forest preserve configuration. :[font = text; inactive; preserveAspect] (3) Step 2 is repeated m times using :[font = input; preserveAspect] NestList[pyro, forestPreserve, m] :[font = text; inactive; preserveAspect; endGroup; endGroup] These pieces of code can be combined into a program. :[font = subsection; inactive; Cclosed; dontNoPageBreakBelow; preserveAspect; startGroup] Program :[font = input; pageBreakBelow; preserveAspect; endGroup; endGroup] smokeyTheBear[n_, p_, f_, m_] := Module[{}, forestPreserve = Table[Random[Integer], {n}]; pyro = (treeGrowIgnite = #/.{0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]}; forestIgnite = TreeGrowIgnite //. {{2, c___, 1} -> {2, c, 2}, {1, c___, 2} -> {2, c, 2}, {a___, 2, 1, b___} -> {a, 2, 2, b}, {a___, 1, 2, b___} -> {a, 2, 2, b}}; forestNew = forestIgnite /. 2 -> 0)&; trees = NestList[pyro, forestPreserve, m] ] :[font = section; inactive; Cclosed; dontNoPageBreakBelow; noKeepOnOnePage; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] Distribution of Forests :[font = text; inactive; preserveAspect] It is useful to see the 'forest through the trees'. To identify the various forests in the forest preserve, a program which labels clusters on a one-dimensional lattice with reflecting boundary conditions can be written consisting of the following steps: :[font = text; inactive; preserveAspect] (1) The lattice sites in the lattice, lat, are updated using :[font = text; inactive; preserveAspect] clusterID[#, RotateLeft[#]]&[lat] :[font = text; inactive; preserveAspect] where :[font = input; preserveAspect] clusterID[1, 0] := i++ clusterID[1, 1] := i clusterID[0, b_] := 0 Attributes[clusterID] = Listable; :[font = text; inactive; preserveAspect] where the two arguments of clusterID are the value of a cell and the value of its right nn. clusterID places an identifying label, i, on the sites in each cluster of contiguous sites. The initial value of i is 1 and increases by 1 for each successive cluster. :[font = text; inactive; preserveAspect] (2) After clusterID has been applied to the lattice, all of the clusters have been correctly labelled with one exception: clusters occurring at the extreme ends of the lattice have been identified as different clusters when they should be treated as a single cluster. To correctly identify the sites in these two clusters as belonging to the same cluster, we can use the anonymous function :[font = input; Cclosed; preserveAspect; startGroup] If[MatchQ[0, First[#] | Last[#]], #, # /. Last[#] -> 1]& :[font = output; output; inactive; preserveAspect; endGroup] If[MatchQ[0, First[#1] | Last[#1]], #1, #1 /. Last[#1] -> 1] & ;[o] If[MatchQ[0, First[#1] | Last[#1]], #1, #1 /. Last[#1] -> 1] & :[font = text; inactive; preserveAspect] These code fragments can be combined into a program. :[font = input; noPageBreak; preserveAspect] clusterLabel1D[lis_List] := Module[{i = 1, clusterID, relabel}, clusterID[1, 0] := i++; clusterID[1, 1] := i; clusterID[0, b_] := 0; Attributes[clusterID] = Listable; result = clusterID[#, RotateLeft[#]]&[lis]; If[MatchQ[0, First[#] | Last[#]], #, # /. Last[#] -> 1]&[result] ] :[font = text; inactive; preserveAspect] The use of the clusterLabel1D function can be be illustrated with a simple forest preserve. :[font = input; Cclosed; preserveAspect; startGroup] config = {1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1}; labelledForests = clusterLabel1D[config] :[font = output; output; inactive; preserveAspect; endGroup] {1, 1, 0, 0, 2, 0, 3, 3, 3, 0, 1} ;[o] {1, 1, 0, 0, 2, 0, 3, 3, 3, 0, 1} :[font = text; inactive; preserveAspect] Having identified the various forests in config, their size distribution can be determined as follows: :[font = text; inactive; preserveAspect] The sizes of the various forests in config is given by :[font = input; Cclosed; preserveAspect; startGroup] forestSizes = Map[Count[labelledForests,#]&, Rest[Union[labelledForests]] ] :[font = output; output; inactive; preserveAspect; endGroup] {3, 1, 3} ;[o] {3, 1, 3} :[font = text; inactive; preserveAspect] The result indicate that the forests in config have 3, 1, and 3 trees, respectively. :[font = text; inactive; preserveAspect] (note: While forestSizes gives the sizes of forests in the order in which they occur in the forest preserve, it does not give the sizes of the deserts lying between the forests. This can be determined using :[font = input; Cclosed; preserveAspect; startGroup] Map[({#, 1})&, config] //. {u___, {v_, r_}, {v_, s_}, w___} -> {u, {v, r + s}, w} /. {{y_, t_}, z___,{y_, u_}}->{{y, t + u}, z} :[font = output; output; inactive; preserveAspect; endGroup] {{1, 3}, {0, 2}, {1, 1}, {0, 1}, {1, 3}, {0, 1}} ;[o] {{1, 3}, {0, 2}, {1, 1}, {0, 1}, {1, 3}, {0, 1}} :[font = text; inactive; preserveAspect] which indicates a sequence consisting of a three-tree forest, a two-site desert, a lone tree, a single empty site, a three-tree forest and an empty site). The distribution of forest sizes is computed using forestSites with the frequency function :[font = input; preserveAspect] frequency[x_List] := Map[{#, Count[x, #]}&, Union[x]] :[font = input; Cclosed; preserveAspect; startGroup] frequency[forestSizes] :[font = output; output; inactive; preserveAspect; endGroup] {{1, 1}, {3, 2}} ;[o] {{1, 1}, {3, 2}} :[font = text; inactive; preserveAspect] which indicates that config contains one lone tree and two 3-tree forests. :[font = text; inactive; preserveAspect; endGroup] A graph of the forest size distribution can be made using ListPlot on the final list. :[font = section; inactive; Cclosed; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] Computer Simulation Projects :[font = text; inactive; preserveAspect] 1. A two-dimensional version of the forest fire algorithm has recently been presented (see Physica A, vol. 204, pp.212-229 (1994)) which appears to show a variety of interesting behaviors. The model is described as follows: ;[s] 2:0,0;3,1;225,-1; 2:1,10,8,Times,1,12,0,0,0;1,11,8,Times,0,12,0,0,0; :[font = text; inactive; preserveAspect] The sites of a square lattice with periodic boundary conditions are either: occupied by trees (ie., have a value 1), burning trees (ie., have a value 2) or are empty (ie., have a value 0). during a time step, all of the sites are updated acccording to the following rules: :[font = text; inactive; preserveAspect] (1) a tree becomes a burning tree with probability (1 - g) if at least one nn is burning. (2) a tree becomes a burning tree with probability f(1 - g) if no nn is burning. (3) a burning tree becomes an empty site. (4) an empty site becomes a tree with probability p. :[font = text; inactive; preserveAspect] where there are three probability parameters: f is the lightning probability, g is the immunity, and p is the tree growth probability. :[font = text; inactive; preserveAspect] It has been found that for this model: (a) when f and g are zero, spiral-shaped fire fronts, reminiscent of the excitable media (see chapter 10), form for small vaules of p. (b) when g = 0, a self-organized critical state (see chapter 8) occurs. (c) when f = 0, a percolation-like phase transition (see chapter 5) takes place at a critical value of g which depends on p. :[font = text; inactive; preserveAspect] The spirals form when empty areas grow at the expense of forests as fire advances while forests grow into empty areas in the absense of fire, so the fire fronts and tree fronts wind around one another. The SOC state results when tree growth occurs more often than lightning and forests burn down faster than trees grow so that fires of all sizes occur. The percolation-like phase transition happens when there is a zero fire density. :[font = text; inactive; preserveAspect] Implement this algorithm in a CA and use it to create snapshots of the forest preserve showing these behaviors. :[font = text; inactive; preserveAspect] 2. The average size of forests in the two-dimensional forest preserve is of interest. To determine this quantity, we can first identify the forests in the preserve and then measue their sizes. ;[s] 2:0,0;2,1;194,-1; 2:1,10,8,Times,1,12,0,0,0;1,11,8,Times,0,12,0,0,0; :[font = text; inactive; preserveAspect] The labeling procedure for the two-dimensional lattice having reflecting boundaries can be explained using a simple lattice as an example. :[font = input; Cclosed; preserveAspect; startGroup] (lat = {{1,0,1,0,0},{0,1,0,1,0},{1,1,0,1,0},{1,0,0,0,1}})//MatrixForm :[font = output; output; inactive; preserveAspect; endGroup] MatrixForm[{{1, 0, 1, 0, 0}, {0, 1, 0, 1, 0}, {1, 1, 0, 1, 0}, {1, 0, 0, 0, 1}}] ;[o] 1 0 1 0 0 0 1 0 1 0 1 1 0 1 0 1 0 0 0 1 :[font = text; inactive; preserveAspect] (1) The lattice sites at the 'northwest' (ie. top left) corners of the clusters are first labelled using :[font = input; preserveAspect] clusterID[RotateRight[lat, {1, 0}], RotateRight[lat, {0, 1}], lat] :[font = text; inactive; preserveAspect] where :[font = input; preserveAspect] clusterCornerID[1, 0, 0] := i++; clusterCornerID[a_ , __] := a; Attributes[clusterCornerID] = Listable :[font = text; inactive; preserveAspect] The three arguments of clusterCornerID are the value of a site, the nn above the site, and the value of the nn to the left of the site. Using clusterCornerID with lat, we get :[font = input; Cclosed; preserveAspect; startGroup] i = 2; clusterCornerID[1,0, 0]:= i++; clusterCornerID[a_, __]:= a; Attributes[clusterCornerID] = Listable; lat = {{1,0,1,0,0},{0,1,0,1,0},{1,1,0,1,0},{1,0,0,0,1}}; (cornerLabels = clusterCornerID[#, RotateRight[#,{1,0}], RotateRight[#,{0,1}] ]&[lat])//MatrixForm :[font = output; output; inactive; preserveAspect; endGroup] MatrixForm[{{1, 0, 2, 0, 0}, {0, 3, 0, 4, 0}, {5, 1, 0, 1, 0}, {1, 0, 0, 0, 6}}] ;[o] 1 0 2 0 0 0 3 0 4 0 5 1 0 1 0 1 0 0 0 6 :[font = text; inactive; preserveAspect] note: In performing the above labeling, the initial value of i was taken to be 2 (and as a result, cluster corner sites were numbered 2, 3, ...) because the value 1 was used to identify non-corner cluster sites. However, i = 1 can be used in the final clusterLabel2D program). :[font = text; inactive; preserveAspect] (2) Having labelled the sites at the corner of clusters, we next label the sites (those with value 1) that lie within clusters and we merge contiguous clusters . Both of these are accomplished with the following rules :[font = input; preserveAspect] reLabel[0, ___] := 0; reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e]; Attributes[reLabel] = Listable; :[font = text; inactive; preserveAspect] where the five arguments of reLabel are: the value of a site, the value of the nn above the site, the value of the nn to the left of the site, the value of the nn beneath the site and the value of the nn to the right of the site.Using FixedPoint to apply reLabel repeatedly to cornerLabels (which is the result of applying clusterCornerID to lat), until the cluster labels no longer change, we get :[font = input; Cclosed; preserveAspect; startGroup] reLabel[0, ___] := 0; reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e]; Attributes[reLabel] = Listable; FixedPoint[reLabel[#, RotateRight[#,{1, 0}], RotateRight[#,{0, 1}], RotateRight[#,{-1,0}], RotateRight[#,{0,-1}]]&, cornerLabels]//MatrixForm :[font = output; output; inactive; preserveAspect; endGroup] MatrixForm[{{6, 0, 2, 0, 0}, {0, 6, 0, 4, 0}, {6, 6, 0, 4, 0}, {6, 0, 0, 0, 6}}] ;[o] 6 0 2 0 0 0 6 0 4 0 6 6 0 4 0 6 0 0 0 6 :[font = text; inactive; preserveAspect] We can now combine these code fragments into a program. :[font = input; preserveAspect] clusterLabel2D[lat_List]:= Module[{i=2}, clusterCornerID[1,0, 0]:= i++; clusterCornerID[a_, __]:= a; Attributes[clusterCornerID] = Listable; cornerLabels = clusterCornerID[#, RotateRight[#,{1,0}], RotateRight[#,{0,1}]]&[lat]; reLabel[0, ___] := 0; reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e]; Attributes[reLabel] = Listable; FixedPoint[reLabel[#, RotateRight[#,{1, 0}], RotateRight[#,{0, 1}], RotateRight[#,{-1,0}], RotateRight[#,{0,-1}]]&, cornerLabels] ] :[font = text; inactive; preserveAspect] If desired, it is straightforward (see the relabelrules2 function in the clusterLabel program in Chapter 5) to relabel the identified clusters so that their numbering is sequential with no gaps. :[font = text; inactive; preserveAspect] While the clusterLabel2D program is specifically designed for use with a lattice system having reflecting boundary conditions, it can be used with a lattice having absorbing boundary conditions by 'decorating' the lattice with rows and columns of 0's (the sandpile model program in Chapter 8 shows how to do this using the absorbBC function ). :[font = text; inactive; preserveAspect; endGroup] Compare the relative speeds of the clusterLabel2D program and the clusterLabel program in chapter 5 as function of the lattice size. :[font = section; inactive; Cclosed; noPageBreak; preserveAspect; fontSize = 14; fontName = "Times"; startGroup] General References :[font = text; inactive; noPageBreak; preserveAspect] Paczuski, Maya and Bak, Per. 1993 "Theory of the one-dimensional forest-fire model", Phys. Rev. E, 48, pp. R3214-R3216. :[font = text; inactive; preserveAspect; endGroup; endGroup] Bak, Per and Paczuski, Maya. 1993 "Why Nature is Complex", Physics World, pp. 396-403. ^*)