Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
1  
Please, provide the code you are working on. –  Sektor Jan 4 at 12:20
1  
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. –  Silvia Jan 4 at 13:07
1  
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? –  Vitaliy Kaurov Jan 4 at 23:16
1  
@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. –  Sophie Jan 5 at 16:02
1  
@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. –  Sophie Jan 5 at 16:07

2 Answers 2

Background

A while back someone told me about a homework problem he was working on. Since I hadn't been in that class I hadn't heard about cellular automatons... So when I was bored and wanted to try the homework problem I ended up implementing everything from scratch. A few weeks later I stumbled upon CellularAutomaton in Vitaliy Kaurov's post here. I hope that my solutions can help you somehow.

The problem

The homework problem is 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 named 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]

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 please know that all I have ever done regarding CA is in this post, but if I can I will try to help you.

share|improve this answer
2  
+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. –  Michael E2 Jan 6 at 6:49

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

share|improve this answer
    
+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. –  Pickett Jan 6 at 10:23
    
@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.... –  Sophie Jan 8 at 21:25
    
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}} .... –  Sophie Jan 8 at 21:25
    
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? –  Sophie Jan 8 at 21:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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