14
$\begingroup$

Im attempting to simulate a simple evacuation model using cellular automata theory.My study aims to determine whether Cellular Automata can effectively represent how people would behave in a fire situation.

I have a grid filled with values that represents the floor plan of a building (see below) and I wish to write rules that will dictate the occupants movement when evacuating. The corridor is supposed to be the safe area although I may have to adjust it slightly so that just the end corridors are the safe area. I am looking to model a fire evacuation, so each cell will contain one person. I was going to assign the value of 0 to cells that contain fire, injured people or debris (obstacles). when designing the simulation by hand I used rules such as if the product of the 2 neighbouring cells was less than 20 the person could move to the next cell, if the product was 20 the cell remained the same and if the product was 0 (when cells meet fire, obstacles etc) then the cell dies and becomes 0 itself. It worked well but obviously I need to now apply a system similar into Mathematica. I have been looking at using 2 options to do this:

  1. use the built in cellular automata function
  2. write my own cell states as functions and use an update function to update the system.

Any suggestions which one would be better to use?

Below is the grid I have set up, the values 500 represent walls, 1 represents the doorways and evacuation route, the other numbers 2-7 represent the distance from the nearest exit.

{
 {500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500},
 {500, 1, 1, 500, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 7, 500, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 500, 5, 5, 5, 5, 5, 500, 1, 1, 500},
 {500, 1, 1, 500, 4, 4, 4, 4, 4, 4, 4, 4, 5, 6, 7, 500, 4, 4, 4, 4, 4, 4, 4, 4, 5, 6, 500, 4, 4, 4, 4, 4, 500, 1, 1, 500},
 {500, 1, 1, 500, 3, 3, 3, 3, 3, 3, 3, 4, 5, 6, 7, 500, 3, 3, 3, 3, 3, 3, 3, 4, 5, 6, 500, 3, 3, 3, 3, 3, 500, 1, 1, 500},
 {500, 1, 1, 500, 2, 2, 2, 2, 2, 2, 3, 4, 5, 6, 7, 500, 2, 2, 2, 2, 2, 2, 3, 4, 5, 6, 500, 3, 2, 2, 2, 3, 500, 1, 1, 500},
 {500, 1, 1, 500, 500, 500, 1, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 1, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 1, 500, 500, 500, 1, 1, 500},
 {500, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 500},
 {500, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 500},
 {500, 1, 1, 500, 500, 500, 500, 500, 500, 500, 500, 500, 1, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 1, 500, 500, 500, 500, 500, 1, 500, 500, 500, 1, 1, 500},
 {500, 1, 1, 500, 7, 6, 5, 4, 3, 2, 2, 2, 2, 2, 2, 500, 6, 5, 4, 3, 2, 2, 2, 2, 2, 2, 500, 3, 2, 2, 2, 3, 500, 1, 1, 500},
 {500, 1, 1, 500, 7, 6, 5, 4, 3, 3, 3, 3, 3, 3, 3, 500, 6, 5, 4, 3, 3, 3, 3, 3, 3, 3, 500, 3, 3, 3, 3, 3, 500, 1, 1, 500},
 {500, 1, 1, 500, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 500, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 500, 4, 4, 4, 4, 4, 500, 1, 1, 500},
 {500, 1, 1, 500, 7, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 500, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 500, 5, 5, 5, 5, 5, 500, 1, 1, 500},
 {500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500}
}

floor evacuation model

$\endgroup$
15
  • 1
    $\begingroup$ Please, provide the code you are working on. $\endgroup$
    – Sektor
    Commented Jan 4, 2014 at 12:20
  • 1
    $\begingroup$ Regarding the rule, have you read the documentation of CellularAutomaton? Especially the Detail -> Possible forms for rule are part. Also, you can have a look at this answer, where I defined a non-trivial custom rule function for CellularAutomaton. $\endgroup$
    – Silvia
    Commented Jan 4, 2014 at 13:07
  • 1
    $\begingroup$ Use @-tag like in @Sophie in comments so we'll get notified. Look through other comments how people do it. What do you expect the program to do and show? Is it a path on the floor? Is it moving people? Is it density of people? $\endgroup$ Commented Jan 4, 2014 at 23:16
  • 1
    $\begingroup$ @VitaliyKaurov the corridor is supposed to be the safe area although I may have to adjust it slightly so that just the end corridors are the safe area. I am looking to model a fire evacuation, so each cell will contain one person. I was going to assign the value of 0 to cells that contain fire, injured people or debris (obstacles). my study is aiming to determine whether Cellular Automata can effectively represent how people would behave in a fire situation. $\endgroup$
    – Sophie
    Commented Jan 5, 2014 at 16:02
  • 1
    $\begingroup$ @VitaliyKaurov when designing the simulation by hand I used rules such as if the product of the 2 neighbouring cells was less than 20 the person could move to the next cell, if the product was 20 the cell remained the same and if the product was 0 (when cells meet fire, obstacles etc) then the cell dies and becomes 0 itself. It worked well but obviously I need to now apply a system similar into Mathematica. $\endgroup$
    – Sophie
    Commented Jan 5, 2014 at 16:07

2 Answers 2

16
$\begingroup$

This is a sample problem that I hope can be helpful; a fire that can spread in four directions (the Von Neumann neighborhood). Only trees can catch fire though. A typical fire looks like this:

forest

Any of the methods below can generate it. There appears to be a name collision between method 1 and 2, so you may need to restart the kernel before you can try the next.

Method #1

(* Generate random forest *)
forest = RandomChoice[{1, 1, 0}, {100, 100}];

(* Set fire to a randomly chosen location *)
forest = ReplacePart[forest, RandomChoice[Position[forest, 1]] -> 3];

(* Kernel. Determines how the fire spreads. *)
ker = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};  (* Can be programmatically generated: Position[CrossMatrix[1],1]-2 *)

(* Find trees in danger *)
trees[forest_] := Position[forest, 1];
dangerZone[forest_] := Intersection[trees[forest], Flatten[ker + ConstantArray[#, Length[ker]] & /@ Position[forest, 3], 1]]

(* Simulates the fire *)
frames = FixedPointList[ReplacePart[#, {Position[#, 3] -> 2, dangerZone[#] -> 3}] &, forest];

(* Visualizes the fire *)
ListAnimate[ ArrayPlot[#, ColorRules -> {0 -> Brown, 1 -> Green, 2 -> Orange, 3 -> Yellow}] & /@ frames]

Method #2

(* Generate random forest *)
forest = RandomChoice[{1, 1, 0}, {100, 100}];

(* Locate trees *)
trees = Position[forest, 1];

(* Set fire to a randomly chosen location *)
forest = ReplacePart[forest, RandomChoice[trees] -> 3];

(* Kernel. Determines how the fire spreads. *)
ker[p_] := Sequence[p + {1, 0}, p + {0, 1}, p + {-1, 0}, p + {0, -1}]; 

(* Find trees in danger *)
dangerZone[frontline_, forest_] := Intersection[trees, DeleteDuplicates[ker /@ frontline]]

(* The brain *)
setFire[frontline_, forest_] := Module[{dz = dangerZone[frontline, forest]},
  If[dz != {},
   trees = Complement[trees, dz];
   setFire[dz, ReplacePart[Sow@forest, {frontline -> 2, dz -> 3}]],
   Sow@forest
   ]
  ]

(* Simulate the fire *)
Block[{$RecursionLimit = 100000},
  {final, {frames}} = Reap@setFire[Position[forest, 3], forest]
  ];

(* Visualizes the fire *)
ListAnimate[ArrayPlot[#, ColorRules -> {0 -> Brown, 1 -> Green, 2 -> Orange, 3 -> Yellow}] & /@ frames]

Method #3

(* Generate random forest *)
forest = ArrayPad[RandomChoice[{1, 1, 0}, {100, 100}], 1];

(* Set fire to a randomly chosen location *)
forest = ReplacePart[forest, RandomChoice[Position[forest, 1]] -> 2];

(* Simulate *)
adv[forest_] := CellularAutomaton[{
    {{_, Except[2], _}, {Except[2], x_, Except[2]}, {_, Except[2], _}} :> x,
    {{_, _, _}, {_, x_, _}, {_, _, _}} :> Switch[x, 1, 2, 2, 2, 0, 0]}, forest];
frames = FixedPointList[adv, forest];

(* Visualize the fire *)
ListAnimate[ArrayPlot[#, ColorRules -> {0 -> Brown, 1 -> Green, 2 -> Orange, 3 -> Yellow}] & /@ frames]

Method #4

New as of the 26th of August 2014 (i.e. not a part of the original post.) This method should be faster than the others.

forest = RandomChoice[{1, 1, 0}, {1000, 1000}];
forest = ReplacePart[forest, RandomChoice[Position[forest, 1]] -> 10];

res = FixedPoint[
   Composition[
    Unitize[First@#] + 9 UnitStep[Last@# - 110] &,
    {#, ListConvolve[{{0, 1, 0}, {1, 100, 1}, {0, 1, 0}}, #, {2, 2}, 0]} &
    ],
   forest
   ];
ArrayPlot[res, ColorRules -> {0 -> Brown, 1 -> Green, 10 -> Orange}]

Conclusions

I hope that in these three different examples you can find something that might help you simulate your model. As for efficiency, a lot of the time is spent visualizing. The computation is much faster, and I tried to measure it, although I am unsure of my results because they seem very unlikely. I would be grateful if someone would confirm them. Remember that these times also seem faster than you anticipate because it takes time for the frontend to display the result, even in list form. For 200x200 pixel grid the methods take this much time:

Method 1: 1.513038 seconds

Method 2: 0.376987 seconds

Method 3: 20-25 seconds

I tried using the third argument of CellularAutomaton to generate only the end state, however it did not yield much improvement in time.

The conclusion is that it appears one can do simulations faster by writing one's own code for it. However CellularAutomaton offers concise syntax and is somewhat easy to use. If you can think of a way to formulate your rules in the format required by CellularAutomaton, therefore, you could start with that and see if speed is a problem or not. In the other methods you are free to write program your rules however you want.

If you don't want to generate the entire simulation at once but just want to watch how it evolves in the beginning, you can save time by visualizing it using Dynamic like this:

Dynamic[ArrayPlot[forest = adv[forest], ColorRules -> {0 -> Brown, 1 -> Green, 2 -> Orange}]]

This might be useful if you are experimenting. Vitaliy uses it in the post I linked to.

If you have any questions you can ask, and I will try to help you.

$\endgroup$
1
  • 2
    $\begingroup$ +1 since I made use of your idea -- by definition this is answer is useful! :P Anyway, nice post. I agree that CA methods are usually not be the fastest. $\endgroup$
    – Michael E2
    Commented Jan 6, 2014 at 6:49
6
$\begingroup$

Here's a way to model @Anon's forest fire with image processing. ImageFilter is not very fast, of which the documentation warns you. However it's acceptable, especially because it makes a nice Dynamic animation. Perhaps it could be adapted - certainly if any of Anon's methods can - to the OP's question. One difference is that in the OP's problem, the building is a constant and there is a second layer representing the people that is dynamic and interacts with the building.


The burn function consists of rules that sets a tree on fire if it is adjacent to a burning tree, and a burning tree burns out.

forest = Colorize[
  Image @ ReplacePart[
    RandomChoice[{1./4, 1./4, 0.}, {100, 100}],
    RandomInteger[{1, 100}, 2] -> 2./4],
  ColorRules -> {0. -> Brown, 1./4 -> Darker@Green, 1./2 -> Red, 3./4 -> Black}];

burn = 
 With[{inert = List @@ Brown // N,
       tree = List @@ Darker@Green // N, 
       fire = List @@ Red // N,
       burned = {0., 0., 0.}},
  Compile[{{m, _Real, 3}},
   If[m[[2, 2]] == inert || m[[2, 2]] == burned, m[[2, 2]],
    If[m[[2, 2]] == fire, burned,
     If[MemberQ[{m[[1, 2]], m[[3, 2]], m[[2, 1]], m[[2, 3]]}, fire],
      fire,
      tree
      ]]],
   RuntimeOptions -> "Speed", CompilationTarget -> "C"]
  ]

Here's a dynamic animation:

Dynamic[forest = ImageFilter[burn, forest, 1, Interleaving -> True]]

Or one can make a movie (just don't use the burned forest - delete the Dynamic output and recompute a new forest!). It took 1.71 seconds for the FixedPointList below.

movie = FixedPointList[ImageFilter[burn, #, 1] &, foo];

Export["foo.gif", movie]

Animation

$\endgroup$
4
  • $\begingroup$ +1. I believe that to adapt it to the OP's problem all one has to do is use the corresponding pixel in the original matrix instead of the color black. It can be done easily for all methods. $\endgroup$
    – C. E.
    Commented Jan 6, 2014 at 10:23
  • $\begingroup$ @Silvia, thank you all for your help I think from all your comments I have successfully computed my cellstates, neighbourhood and rules however when finally attempting to run the cellular automaton function I am getting this note.... $\endgroup$
    – Sophie
    Commented Jan 8, 2014 at 21:25
  • $\begingroup$ CellularAutomaton::initn: The initial condition specification {{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},{1,0,0,1,5,2,2,0,2,2,0,0,2,0,2,1,2,0,0,0,2,0,2,2,2,0,1,2,0,0,2,0,1,0,0,1},{1,0,0,1,0,2,0,2,0,0,2,2,0,2,0,1,0,2,2,2,0,2,2,0,2,0,1,2,0,0,2,0,1,0,0,1},<<8>>,{1,0,0,1,2,2,2,2,2,2,2,2,2,2,0,1,2,2,2,2,0,2,2,2,2,0,1,0,2,0,2,0,1,0,0,1},{1,0,0,1,2,2,2,2,0,2,2,0,2,0,0,1,2,0,0,2,2,2,2,2,2,0,1,0,2,0,2,0,1,0,0,1},{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}} .... $\endgroup$
    – Sophie
    Commented Jan 8, 2014 at 21:25
  • $\begingroup$ should be of the form aspec, {aspec, bspec}, or {{{aspec1, off1}, {aspec2, off2},... {aspecn, offn}}, bspec} (n > 0). Each aspec must be a non-empty rank 1 array whose elements at level 1 are integers i in the range 0 <= i <= 2. >> I am going over the cellular automata intial conditions details now but Im worried it wont accept my array which is essential that it remains as it represents the building layout. Any advice? $\endgroup$
    – Sophie
    Commented Jan 8, 2014 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.