Techniques to understand computer simulations: Markov chain analysis. 
Figure 11

Luis R. Izquierdo, Segismundo S. Izquierdo, José M. Galán & José I. Santos

Program to create figures showing the transient distribution of the location of a single walker in a CoolWorld model where the temperature profile and the distribution of houses are as described in section 5.2. The height of each patch denotes the relative probability that the walker is on the patch in any given time-step.

Licence

Copyright (C) 2008 Luis R. Izquierdo, Segismundo S. Izquierdo, José M. Galán & José I. Santos

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110 - 1301, USA.

Contact information :
Luis R. Izquierdo
University of Burgos, Spain.
e - mail : luis@izquierdo.name

Setting the Environment

SetDirectory[NotebookDirectory[]];

Parameters

probRandomMove=0.5;
probLeavingHome=0.01;
xDimension=33;
xMin=-16;
yDimension=33;
yMin=-16;

map=Import["map.png"];

(* Timesteps for which we want a figure showing the transient distribution. Initial conditions are always shown. *)
reportingList=Join[Table[i,{i,1,50}],Table[i,{i,55,195,5}],Table[i,{i,200,300,10}],Table[i,{i,350,500,50}],Table[i,{i,600,1000,100}],Table[i,{i,2000,5000,1000}]];

Function to create a figure from a probability distribution over the state space

numberOfDigits=Length[IntegerDigits[Last[reportingList]]];
points=map[[1,1]]/255.;
ticks=Table[{0.05 k,ToString[NumberForm[0.05 k,{2,2}]]},{k,0,20}];

CreateFigure[stateHistogram_List,timestep_]:=Module[{stateMatrix=N[Partition[stateHistogram,33]],my3DGraph},
stateMatrix=ArrayFlatten[MapIndexed[ConstantArray[#,{13,13}]&,stateMatrix,{2}]];

my3DGraph=ListPlot3D[stateMatrix,Mesh->None,VertexColors->{points},PlotRange->{0,0.25},ViewPoint->{1,-1,0.5},Axes->{False,False,True},Boxed->False,Ticks->{None,None,ticks},
Epilog->Inset[Framed[PaddedForm[timestep,numberOfDigits],Background->LightYellow],{1,0.7},BaseStyle->"Title",Alignment->Right]];

Export[StringJoin[ToString[timestep],".gif"],my3DGraph,ImageSize->800];
];

Creation of the transition matrix

TRANSLATION: state < - > coordinates

offset=-(xMin+xDimension*yMin)+1;
state[{x_,y_}]:=x+xDimension*y+offset;
coordinates[s_]:={Mod[s-1,xDimension]+xMin,Floor[(s-1)/xDimension]+yMin}

TEMPERATURE-DRIVEN MOVE

HottestNeighbour[{x_,y_}]:={x-Sign[x],y-Sign[y]}
HottestState[s_]:=state[HottestNeighbour[coordinates[s]]]

(*temperatureTM, for temperature Transition Matrix, is a transition matrix that characterises the deterministic movement from each patch to its hottest neighbour.*)
temperatureTM=SparseArray[Map[Rule[{#,HottestState[#]},1]&,Range[xDimension*yDimension]],{xDimension*yDimension,xDimension*yDimension}];

RANDOM MOVE

Neighbours[{x_,y_}]:=DeleteCases[Flatten[Table[{i,j},{i,x-1,x+1},{j,y-1,y+1}],1],{α_,β_}/;Or[α<xMin,α≥xDimension+xMin,β<yMin,β≥yDimension+yMin,And[α==x,β==y]]];
Neighbours[{x1_,y1_},{x2_,y2_}]:=If[Abs[x1-x2]≤1&&Abs[y1-y2]≤1&&Or[x1≠x2,y1≠y2],1,0];
Neighbours[s1_,s2_]:=Neighbours[coordinates[s1],coordinates[s2]]
Neighbours[s_]:=Map[state,Neighbours[coordinates[s]]]
RandomMoveRules[s_]:=Module[{n=Neighbours[s],l},
l=Length[n];
Map[Rule[{s,#},1/l]&,n]
]

(*randMoveTM, for random Move Transition Matrix, is a transition matrix that characterises the random movement from one patch to one of its neighbours.*)
randMoveTM=SparseArray[Flatten[Map[RandomMoveRules,Range[xDimension*yDimension]]],{xDimension*yDimension,xDimension*yDimension}];

IN-HOUSE MOVE

patchesWithHouse=Cases[Flatten[Table[{i,j},{i,-5,5},{j,-5,5}],1],v_/;5≤Norm[v]<6];

(*rdMoveRulesForHouses is a set of rules that characterise the random movement from each patch with a house to one of its neighbours.*)
rdMoveRulesForHouses=Flatten[Map[RandomMoveRules[state[#]]&,patchesWithHouse]];

(*houseMoveTM, for in-house Move Transition Matrix, is a transition matrix that characterises the transitions that may occur from one state with a house.*)
houseMoveTM=probLeavingHome*SparseArray[rdMoveRulesForHouses,{xDimension*yDimension,xDimension*yDimension}]+(1-probLeavingHome)*SparseArray[Map[Rule[{state[#],state[#]},1]&,patchesWithHouse],{xDimension*yDimension,xDimension*yDimension}];

FINAL TRANSITION MATRIX

transitionMatrix=SparseArray[Join[Most[ArrayRules[houseMoveTM]],ArrayRules[probRandomMove*randMoveTM + (1-probRandomMove)*temperatureTM]],{xDimension*yDimension,xDimension*yDimension}];

Transient Distributions

numberOfSteps=Last[reportingList];
reportCounter=1;

d1=StringJoin[NotebookDirectory[],"fig11"];
CreateDirectory[d1];
SetDirectory[d1];

(* initial state *)
timestep=0;
x=UnitVector[xDimension*yDimension,1];
CreateFigure[x,timestep++];

(* transient distributions *)
While[timestep≤numberOfSteps,

x=(x.transitionMatrix);

If[(timestep==reportingList[[reportCounter]]),
reportCounter++;
CreateFigure[x,timestep];
];

timestep++;
]

Asymptotic Distribution

x=x.MatrixPower[N[transitionMatrix],100000];
CreateFigure[x,∞];