## Thursday, June 14, 2012

### Werewolf attack: Spatial Multi-agent simulation - basic artificial intelligence

* This is a simulation of a werewolf attack using Stata.

* Initially there is X number of people on a map and a few werewolves.

* Werewolves head toward whatever unprotected humans are on the map.

* If they encounter a human they kill that human and go towards they next unprotected human.

* Humans do not know about the werewolves until they become within sensory range - specified by the user.

* Once they become within sensory range humans immediately seek the closest shelter.

* Shelters are distributed throughout the map.

* Once a human is in a shelter the human is safe.

***************************************************************
* Simulation Parameters Begin
***************************************************************

clear

* Set the number of humans
local num_humans=1500

local human_vision=65

* Human run speed
local human_speed = 1

* Set number of shelters
local num_shelters = 3

* Specify the initial number of werewolves
local num_werewolves = 8

* Werewolf run speed (if humans are faster than werewolf then it is only the unlucky human that will ever be caught)
local werewolf_speed = 5

* Specify the dimensions of the world grid
local xrange=400
local yrange=400

* Graph interval every x rounds the program graphs what happens
local graph_int = 20

* Tell stata to not graph as the simulation progresses
local draw=1

***************************************************************
* Simulation Parameters End
***************************************************************

* This tells stata to either draw the graphs as it runs the code or wait till the end to draw the graphs.
if `draw'==0 local nodraw nodraw

set obs `num_humans'

gen hmn_id=_n
label var hmn_id "Human ID"

* To begin with all types are human
gen ty = 1
lab var ty "Observation type 1=human, 2=werewolf, 3=shelter"
label define obj_type 1 "Human" 2 "Werewolf" 3 "Shelter" 4 ///
"Deceased" 5 "Spooked Humans" 6 "Safe"
label values ty obj_type

* Then we will add the werewolves

set obs `=`num_humans'+`num_werewolves''

replace ty = 2 if ty==.

* Put the werewolves on the top of the list
gsort -ty
gen were_id = _n if ty==2
label var were_id "Werewolf ID"

* Next the shelters
set obs `=`num_humans'+`num_werewolves'+`num_shelters''

replace ty = 3 if ty==.

tab ty
gsort -ty
gen shelter_id = _n if ty==3
label var shelter_id "Shelter ID"

* Looks pretty good so far.

* People shelters and werewolves are uniformly throughout the data.
gen x=runiform()*`xrange'
gen y=runiform()*`yrange'

* Though werewolves will enter the maps from the edges
gen xy_edge = rbinomial(1,.5)

replace x=round(x, `xrange') if xy_edge==0 & ty == 2
replace y=round(y, `yrange') if xy_edge==1 & ty == 2
* This makes it so the werewolves start on one of the sides of the map.

* Generate
sum x
replace x=rnormal(r(mean),r(sd)*5/8) if ty == 3
sum y
replace y=rnormal(r(mean),r(sd)*5/8) if ty == 3

two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
(scatter y x if ty==2, color(cranberry) )               ///
, legend(label (1 "Humans")                             ///
label (2 "Werewolves")                                   ///
label (3 "Shelters")                                    ///
rows(1))  `nodraw' title(Full Moon - Round 0)           ///
plotregion(fcolor(gs2))

***************************************************************
* First we will calculate which direction and where the humans will run.
* This direction is constant since humans do not try to avoid werewolves, they just run for shelter.
***************************************************************

* We will loops through shelters and calculate distance for each human.

global shelter_distance_list

qui forv i=1(1)`num_shelters' {
* The first x number of observations are the werewolves.
sort shelter_id
local shelter_x = x[`i']
local shelter_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey

* Now let's calculate the distance of the spooked humans from the shelters.
cap drop dist_hm_shltr`i'
gen dist_hm_shltr`i'=((x-`shelter_x')^2 + (y-`shelter_y')^2)^.5 if ty==1

global shelter_distance_list \$shelter_distance_list , dist_hm_shltr`i'

}

di "\$shelter_distance_list"
* We can see in out list of shelters there is unwanted comma.

* For now we will just fill in an meaningless value to prevent error.

* Now we figure out which shelter is the closest for all of the humans collectively
gen closest_shelter = .
gen closest_shelter_x = .
gen closest_shelter_y = .
gen closest_shelter_dist = .

qui forv i=1(1)`num_shelters' {
replace closest_shelter = `i' if  dist_hm_shltr`i' == min(99999 \$shelter_distance_list)
replace closest_shelter_dist = dist_hm_shltr`i' if  closest_shelter == `i'

local shelter_x = x[`i']
local shelter_y = y[`i']

replace closest_shelter_x = `shelter_x' if closest_shelter == `i'
replace closest_shelter_y = `shelter_y' if closest_shelter == `i'

}

tab closest_shelter

* Now we will calculate the ideal trajectory for each human given that human is running for shelter.
* We know a, b, c [see graph] and the next section for algebra.
gen a = y-closest_shelter_y
gen b = x-closest_shelter_x

gen a_prime=1/(1+b^2/a^2)^.5

* b'^2 = 1-a'^2
* b' = (1-a'^2)^.5
gen b_prime = (1-(a_prime)^2)^.5

gen x2 = x - b_prime*sign(b) if ty == 1
gen y2 = y - a_prime*sign(a) if ty == 1

replace x2 = x2 - b_prime*sign(b) if ty == 1
replace y2 = y2 - a_prime*sign(a) if ty == 1

two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y2 x2 if ty==1, mcolor(sandb) msymbol(smcircle)) ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
, legend(label (1 "Starting Place")                     ///
label (2 "One step")                                    ///
label (3 "Shelters")                                    ///
rows(1)) `nodraw' title(Sample humans running for shelter) ///
legend(color(white) region(fcolor(black)))              ///
plotregion(fcolor(gs2))

drop x2 y2

***************************************************************
* Werewolves search for human targets
***************************************************************

* I loop through each of the werewolves.

gen x2=.
gen y2=.

* Each werewolf will survey all of the humans and find the one who is the closest.
qui forv i=1(1)`num_werewolves' {
* The first x number of observations are the werewolves.
sort were_id
local were_x = x[`i']
local were_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey

* Now let's calculate the distance of that werewolf from the humans.
cap drop dist_hm_wr
gen dist_hm_wr`i'=((x-`were_x')^2 + (y-`were_y')^2)^.5 if ty==1 | ty==5
* We take the square but for finding targets it is unnecessary since the ranks of the different distances remains constant.

replace ty=5 if dist_hm_wr`i'  <  `human_vision'

sort dist_hm_wr
* The top of this list is the target human.

* If a human is less than `werewolf_speed' distance from the were and it is the closest to that were then the human is dead.
replace ty = 4 if dist_hm_wr < `werewolf_speed' &  _n==1

* Now we need to calculate the attack vector.

* Ie, what is the linear projection that gets the werewolf from its current location to the human fastest (in terms of x and y movement)?

* We can use some algebra and a little trig to find the proper vector.

* We know a, b, c
local a = `were_y'-y[1]
local b = `were_x'-x[1]

* And c' which is equal to the speed of werewolves
local c_prime=`werewolf_speed'

* Triangle abc is similar to triangle a'b'c' so the ratios b/a=b'/a'
* b'=b(a'/a)
* (1) b'^2=b^2 * a'^2 / a^2
* (2) Pathagoras (500BC): a'^2 + b'^2 = c'^2
* (3) instert (1) into (2): a'^2 + b^2 * a'^2 / a^2 = c'^2
* a'^2(1 + b^2 / a^2) = c'^2
* a'^2 = c'^2/(1 + b^2 / a^2)
* a'=c'/(1 + b^2 / a^2)^.5

local a_prime=`c_prime'/(1 + (`b')^2 / (`a')^2)^.5

* b'^2 = c'^2-a'^2
* b' = (c'^2-a'^2)^.5
local b_prime = ((`c_prime')^2-(`a_prime')^2)^.5

* It is important to note the direction that the werewolf needs to travel as well.

/*
replace x = x + `b_prime'*sign(`b') if were_id == `i'
replace y = y + `a_prime'*sign(`a') if were_id == `i'
*/
replace x2 = x - `b_prime'*sign(`b') if were_id == `i'
replace y2 = y - `a_prime'*sign(`a') if were_id == `i'

noi di "Werewolf `i' moves"
}

two (scatter y x if ty==1 , mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==5 , color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry))                ///
(scatter y2 x2 if ty==2, color(red))                    ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
(scatter y x if ty==4, msymbol(smplus) color(gs16))     ///
, legend(label (1 "Unconcerned Humans")                 ///
label (2 "Spooked Humans")                              ///
label (3 "Werewolves")                                   ///
label (4 "Werewolf's 1st Move")                         ///
label (5 "Shelters")                                    ///
label (6 "Deceased")                                    ///
rows(2)) `nodraw' title(Full Moon - Round 0)            ///
legend(color(white) region(fcolor(black)))              ///
plotregion(fcolor(gs2))

***************************************************************
* Spooked Humans Seek Out Closest Shelter
***************************************************************

* Move spooked humans towards shelters.
replace closest_shelter_dist = closest_shelter_dist-`human_speed'

replace ty=6 if closest_shelter_dist < 0
* Once a human gets to a shelter, that human is considered safe.

replace x2 = x - b_prime*sign(b)*`human_speed' if ty == 5
replace y2 = y - a_prime*sign(a)*`human_speed' if ty == 5

two (scatter y x if ty==1 | ty==5, mcolor(gs12) msymbol(smcircle)) ///
(scatter y2 x2 if ty==5 , color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry))                ///
(scatter y2 x2 if ty==2, color(red))                    ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
(scatter y x if ty==4, msymbol(smplus) color(gs16))     ///
, legend(label (1 "Unconcerned Humans")                 ///
label (2 "Spooked Humans")                              ///
label (3 "Werewolves")                                   ///
label (4 "Werewolf's 1st Move")                         ///
label (5 "Shelters")                                    ///
label (6 "Deceased")                                    ///
rows(2)) `nodraw' title(Full Moon - Round 0)            ///
legend(color(white) region(fcolor(black)))              ///
plotregion(fcolor(gs2))

**** LOOKS like things are working properly
* Let's make the position changes to x and y permanent

replace x = x2 if ty == 5 | ty == 2
replace y = y2 if ty == 5 | ty == 2

drop x2 y2

***************************************************************
* Let's allow for some summary statistics

gen obs_num = _n

gen unaware_humans = `num_humans' if _n == 1
gen spooked_humans = 0            if _n == 1
gen deceased_humans = 0           if _n == 1
gen safe_human = 0                if _n == 1

***************************************************************
* Now we will begin the system loop
***************************************************************
global graph_list

cap
local ii=0
while _rc==0 {

local ii=`ii'+1

noi di "Round `=`ii'+1'"
***************************************************************
* Let's record some summary statistics

qui sum ty if ty==1
replace unaware_humans = r(N) if obs_num == `ii'+1

qui sum ty if ty==5
replace spooked_humans = r(N) if obs_num == `ii'+1

qui sum ty if ty==4
replace deceased_humans = r(N) if obs_num == `ii'+1

qui sum ty if ty==6
replace safe_human = r(N) if obs_num == `ii'+1

***************************************************************
* First the werewolf moves

* Each werewolf will survey all of the humans and find the one who is the closest.
sum ty if ty==1 | ty==5
if r(N)==0 cap end_loop
if r(N)>0 qui forv i=1(1)`num_werewolves' {
* The first x number of observations are the werewolves.
sort were_id
local were_x = x[`i']
local were_y = y[`i']
* were_x and were_y hold the position of the current werewolf looking for its prey

* Now let's calculate the distance of that werewolf from the humans.
cap drop dist_hm_wr
gen dist_hm_wr`i'=((x-`were_x')^2 + (y-`were_y')^2)^.5 if ty==1 | ty==5
* We take the square but for finding targets it is unnecessary since the ranks of the different distances remains constant.

replace ty=5 if dist_hm_wr`i' < `human_vision'

sort dist_hm_wr
* The top of this list is the target human.

* If a human is less than `werewolf_speed' distance from the were and it is the closest to that were then the human is dead.
replace ty = 4 if dist_hm_wr < `werewolf_speed' & _n==1

* Now we need to calculate the attack vector.

* Ie, what is the linear projection that gets the werewolf from its current location to the human fastest (in terms of x and y movement)?

* We can use some algebra and a little trig to find the proper vector.

* We know a, b, c
local a = `were_y'-y[1]
local b = `were_x'-x[1]

* And c' which is equal to the speed of werewolves.
local c_prime=`werewolf_speed'

* Triangle abc is similar to triangle a'b'c' so the ratios b/a=b'/a'
* b'=b(a'/a)
* (1) b'^2=b^2 * a'^2 / a^2
* (2) Pathagoras (500BC): a'^2 + b'^2 = c'^2
* (3) insert (1) into (2): a'^2 + b^2 * a'^2 / a^2 = c'^2
* a'^2(1 + b^2 / a^2) = c'^2
* a'^2 = c'^2/(1 + b^2 / a^2)
* a'=c'/(1 + b^2 / a^2)^.5

local a_prime=`c_prime'/(1 + (`b')^2 / (`a')^2)^.5

* b'^2 = c'^2-a'^2
* b' = (c'^2-a'^2)^.5
local b_prime = ((`c_prime')^2-(`a_prime')^2)^.5

replace x = x - `b_prime'*sign(`b') if were_id == `i'
replace y = y - `a_prime'*sign(`a') if were_id == `i'
}

***************************************************************
* Then the humans move

* Move spooked humans towards shelters.
replace closest_shelter_dist = closest_shelter_dist-`human_speed' if ty == 5

replace ty=6 if closest_shelter_dist < 0 & ty == 5
* Once a human gets to a shelter, that human is considered safe.

replace x = x - b_prime*sign(b)*`human_speed' if ty == 5
replace y = y - a_prime*sign(a)*`human_speed' if ty == 5

***************************************************************
* Then the humans move
if `ii' / `graph_int'==int(`ii' / `graph_int') {

two (scatter y x if ty==1, mcolor(gs12) msymbol(smcircle)) ///
(scatter y x if ty==5, color(sandb) msymbol(smcircle)) ///
(scatter y x if ty==2, color(cranberry))                ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
(scatter y x if ty==4, msymbol(smplus) color(gs16))     ///
, legend(label (1 "Unconcerned Humans")                 ///
label (2 "Spooked Humans")                              ///
label (3 "Werewolves")                                   ///
label (4 "Shelters")                                    ///
label (5 "Deceased")                                    ///
rows(2)) `nodraw' title(Full Moon - Round `ii')          ///
legend(color(white) region(fcolor(black)))              ///
plotregion(fcolor(gs2))

two (scatter y x if ty==1, mcolor(gs12) msymbol(smcircle))  ///
(scatter y x if ty==5, color(sandb) msymbol(smcircle))  ///
(scatter y x if ty==2, color(cranberry))                ///
(scatter y x if ty==3, msymbol(square) color(yellow))   ///
(scatter y x if ty==4, msymbol(smplus) color(gs16))     ///
, legend(off) nodraw name(round`ii', replace)           ///
plotregion(fcolor(gs2))

global graph_list \$graph_list round`ii'

}
}

***************************************************************
* Now we will end the system loop
***************************************************************

sort obs_num

gen round=obs_num if unaware_humans!=.

two (line unaware_humans round, sort)                       ///
(line spooked_humans round, sort)                       ///
(line deceased_humans round, sort)                      ///
(line safe_human round, sort)                           ///
, legend(label (1 "Unconcerned Humans")                 ///
label (2 "Spooked Humans")                              ///
label (3 "Deceased Humans")                             ///
label (4 "Safe Humans"))

* Different specifications at the beginning of the model yield different results. (werewolves always win though)
sleep 5000

graph combine \$graph_list

* (c) Francis Smart 2012.