View source for OR/16/022 Appendix 1 - NetLogo associated scripts

__NOTOC__
{{OR/16/022}}
'''NetLogo Script for the Northampton model'''

''NetLogo Script''

extensions [gis]

globals [elevation-dataset]

patches-own [elevation]

turtles-own [exposure group]

to setup

let xHm [-7 6 -5 3]

let yHm [6 5 -5 -6]

set-default-shape turtles "person"

set        elevation-dataset        gis:load-dataset             "C:/Users/mrca/My Documents/Nthpton_As.asc"

gis:paint elevation-dataset 0

gis:apply-raster elevation-dataset elevation

let min-elevation gis:minimum-of elevation-dataset

let max-elevation gis:maximum-of elevation-dataset

ask patches

[ if (elevation <= 0) or (elevation >= 0)

[set pcolor scale-color black elevation min-elevation max-elevation]

]

ask patches

[ if pcolor = black

[set elevation 0]

]

create-turtles numturtles [

set exposure 0

set xcor -7

set ycor 6

set group 1

set color red

]

create-turtles numturtles [

set exposure 0

set xcor 6

set ycor 5

set group 2

set color green

]

create-turtles numturtles [

set exposure 0

set xcor -5

set ycor -5

set group 3

set color blue

]

create-turtles numturtles [

set exposure 0

set xcor 3

set ycor -6

set group 4

set color yellow

]

end to go

ask turtles with [group = 1][

ifelse distancexy -7 6 > dist-to-home

[facexy -7 6]

[rt random 360]               ; set random heading

forward 1           ; advance one step

]

ask turtles with [group = 2][

ifelse distancexy 6 5 > dist-to-home

[facexy 6 5]

[rt random 360]               ; set random heading

forward 1           ; advance one step

]

ask turtles with [group = 3][

ifelse distancexy -5 -5 > dist-to-home

[facexy -5 -5]

[rt random 360]               ; set random heading

forward 1           ; advance one step

]

ask turtles with [group = 4][

ifelse distancexy 3 -6 > dist-to-home

[facexy 3 -6]

[rt random 360]               ; set random heading

forward 1           ; advance one step

]

ask patches [

let my-turtles turtles-here

if any? my-turtles [

let lvls elevation / 100

; print lvls

ask my-turtles

[set exposure (exposure + lvls)]

]

]

End.


'''R Script for the NetLogo Northampton Model'''

library(RNetLogo)

<nowiki>#</nowiki> this part starts up the NetLogo language

nl.path <-"C:/Program Files/NetLogo 5.3/app"

NLStart(nl.path,gui=TRUE, nl.obj=NULL, is3d=FALSE)

<nowiki>#</nowiki> now get the file path for the model we ant to run

model.path<-file.choose()

NLLoadModel(model.path)

<nowiki>#</nowiki> This part sets up the model setting the number of turtles at each site

<nowiki>#</nowiki> and the maximum distance from home the agents are allowed to move

NLCommand("clear-all")

NLCommand("set numturtles 100")

NLCommand("set dist-to-home 3")

NLCommand("setup")

<nowiki>#</nowiki> Now we need to set up some dataframes to store the arsenic exposure

<nowiki>#</nowiki> for each agent in each group

df1<-data.frame(nms=as.character(1:100))

df2<-data.frame(nms=as.character(1:100))

df3<-data.frame(nms=as.character(1:100))

df4<-data.frame(nms=as.character(1:100))

<nowiki>#</nowiki> this loop runs the NetLogo Program over 100 cycles collecting the exposure

<nowiki>#</nowiki> from each agent in each group at each cycle

for(i in 1:100){

NLDoCommand(1, "go")

Exp1<- NLGetAgentSet("exposure","turtles with [group = 1]")

df1<-cbind(df1,Exp1)

Exp2<- NLGetAgentSet("exposure","turtles with [ group = 2 ]")

df2<-cbind(df2,Exp2)

Exp3<- NLGetAgentSet("exposure","turtles with [ group = 3 ]")

df3<-cbind(df3,Exp3)

Exp4<- NLGetAgentSet("exposure","turtles with [ group = 4 ]")

df4<-cbind(df4,Exp4)

}

<nowiki>#</nowiki> now the results can be plotted out

<nowiki>###################################</nowiki>

<nowiki>#</nowiki> This plots the cumulative expusure for each agent in each group

<nowiki>#</nowiki> along with an average value for each group

par(mfrow=c(2,2))

avg1<-sapply(df1[,2:101],FUN=mean)

plot(1:100,avg1,type="l", ylim=c(0,3.5),lty=2,lwd=4,

xlab="time",ylab="Exposure", main="Group 1")

for( i in 1:100){

lines(1:100,df1[i,2:101],col=83+i)

}

lines(1:100,avg1,lty=2,lwd=3)

avg2<-sapply(df2[,2:101],FUN=mean)

plot(1:100,avg2,type="l", ylim=c(0,3.5),lty=2,lwd=2,

xlab="time",ylab="Exposure", main="Group 2")

for( i in 1:100){

lines(1:100,df2[i,2:101],col=83+i)

}

lines(1:100,avg2,lty=2,lwd=3)

avg3<-sapply(df3[,2:101],FUN=mean)

plot(1:100,avg3,type="l", ylim=c(0,3.5),lty=2,lwd=2,

xlab="time",ylab="Exposure", main="Group 3")

for( i in 1:100){

lines(1:100,df3[i,2:101],col=83+i)

}

lines(1:100,avg3,lty=2,lwd=3)

avg4<-sapply(df4[,2:101],FUN=mean)

plot(1:100,avg4,type="l", ylim=c(0,3.5),lty=2,lwd=2,

xlab="time",ylab="Exposure", main="Group 4")

for( i in 1:100){

lines(1:100,df4[i,2:101],col=83+i)

}

lines(1:100,avg4,lty=2,lwd=3)

<nowiki>###################################################</nowiki>

<nowiki>#</nowiki> This does a box and whisker plot for the final exposures over

<nowiki>#</nowiki> 100 cycles par(mfrow=c(1,1))

alldat<-cbind(df1[,101],df2[,101],df3[,101],df4[,101])

boxplot(alldat, xlab="Group",ylab="As Exposure")

<nowiki>###################################################</nowiki>

<nowiki>#</nowiki> This does a probability density function of the final exposure in each

<nowiki>#</nowiki> group

plot(density(alldat[,1]),main="Exposure PDF",

xlab="Cumulative As Exposure",col="green",

ylim=c(0,2),xlim=c(1.5,4.5),lwd=2)

lines(density(alldat[,2]),col="red",lwd=2)

lines(density(alldat[,3]),col="blue",lwd=2)

lines(density(alldat[,4]),col="purple",lwd=2)

legend("topleft",c("Group1","Group2","Group3","Group4"),

col=c("green","red","blue","purple"),bty="n",lty=1,lwd=2)

<nowiki>#################################################</nowiki>

<nowiki>#</nowiki> this quits the NetLogo program

NLQuit()

.
[[Category: OR/16/022 Agent Based Modelling: Initial assessment for use on soil bioaccessibility | 07]]