# Function to create random matrix, then a simulation. *EDIT* Create two functions

How can I write a function with parameters b, i, and j, which would return a i x j array of 0s, 1s, and 2s, but the array needs to be populated randomly with 0s and 1s which have respective probabilities 1/(1+b) and b/(1+b) and then one square selected at random should be converted to a 2? What I have done so far is this,

initialState[b_, i_, _] := RandomChoice[{1/(1 + b), b/(1 + b)} -> {0,1}, {i,j}]


but I don't know how I can select a specific square at random to convert it to 2. How can I go about doing this? And then, how can this function be used to simulate a 30x30 fire, which would start at the array point that is equal to 2, with values of b equal to 0.6 up to 1.8, in increments of 0.2?

EDIT: As a further development of this, how could I write the following two functions: a function that generates the number of trees that have survived after the fire, which run from initial state mat and a second function that generates which trees (represented by 1s) don't have any trees, burning or not, in their von Neumann neighborhood? Then for a 30x30 matrix and varying values of parameter b, how can I calculate a number of iterations for the two functions and then get the average of each, for the different values of b?

Update: Animation of fire-spreading:

ClearAll[initState, vNNeighbors, step, iterationList]

initState[b_, i_, j_, pos_: Automatic] := ReplacePart[
RandomChoice[{1/(1 + b), b/(1 + b)} -> {0, 1}, {i, j}],
(pos /. Automatic -> RandomChoice[Tuples[{Range@i, Range@j}]]) -> 2]

vNNeighbors[dim_: {30, 30}] := AdjacencyList[NearestNeighborGraph@Tuples@Range@dim, #]&

step = MapAt[Min[2, 2 #] &, #, vNNeighbors[][Position[#, 2]]] &;

iterationList[nmax_: Automatic][b_, i_, j_, pos_: Automatic] :=
NestList[step, initState[b, i, j, pos], nmax /. Automatic -> Max[i, j]]


Example:

SeedRandom[1]
startpos = {15, 5};
{i, j} = {30, 30};

Manipulate[ListAnimate[MatrixPlot[#, Mesh -> All, Frame -> False] & /@
iterationList[][b, i, j, startpos], Paneled -> False],
{{b, .6}, .6, 1.6, .2}]


Alternatively, you can use FixedPointList to generate the list of arrrays:

ClearAll[FPList]
FPList[b_, i_, j_, pos_: Automatic] := FixedPointList[step, initState[b, i, j, pos]]

SeedRandom[1]
Manipulate[ListAnimate[MatrixPlot[#, Mesh -> All, Frame -> False] & /@
FPList[b, i, j, startpos], Paneled -> False],
{{b, .6}, .6, 1.6, .2}]


ClearAll[initState]
initState[b_, i_, j_, pos_: Automatic] := ReplacePart[
RandomChoice[{1/(1 + b), b/(1 + b)} -> {0, 1}, {i, j}],
(pos /. Automatic -> RandomChoice[Tuples[{Range@i, Range@j}]]) -> 2]


Default position of 2 is random:

SeedRandom[1]
MatrixPlot[initState[.6, 30, 30], Mesh -> All]


Position of 2 is {10, 15}

SeedRandom[1]
MatrixPlot[initState[.6, 30, 30, {10, 15}], Mesh -> All]


Animations:

Position of 2 changes randomly in every iteration:

SeedRandom[1]
frames1 = Table[MatrixPlot[initState[b, 30, 30], Mesh -> All,
Frame -> False], {b, Range[.6, 1.8, .2]}];

Export["animation1.gif", frames1]


Position of 2 is remains fixed at a random value over iterations:

SeedRandom[1]
frames2 = With[{pos2 = RandomChoice[Tuples[{Range@30, Range@30}]]},
Table[MatrixPlot[initState[b, 30, 30, pos2], Mesh -> All,
Frame -> False], {b, Range[.6, 1.8, .2]}]];

Export["animation2.gif", frames2]


• thank you, this is very helpful. Using ListAnimate, how could I make it so as the 1s in the von Neumann neighborhood of the 2, also turn into 2 as the animation progresses? For the same b parameter values. The animation would simulate a fire spreading, or the 1s turning into 2s
– jbl
Commented Feb 5, 2021 at 20:16
• @jbl, please see the update.
– kglr
Commented Feb 5, 2021 at 22:20

Well, you can sample two random numbers and use them as indices of your array:

initialState[b_, i_, j_] :=
Module[{i2 = RandomInteger[{1, i}], j2 = RandomInteger[{1, j}],
randArray =
RandomChoice[{1/(1 + b), b/(1 + b)} -> {0, 1}, {i, j}]},
randArray[[i2, j2]] = 2; randArray]


Unfortunately, I do not undestand the rest of your question.

• apologies for not phrasing it properly. Basically, I need to create a ListAnimate plot which simulates a forest fire of a 30 x 30 array, in which the fire starts where 2 is at in the array. This then spreads to its von Neumann neighbors that are equal to 1. Essentially, the 1s turn into 2s as the fire spreads. Sorry if this still not very clear
– jbl
Commented Feb 5, 2021 at 19:48